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Chapter  1 
INTRODUCTION 

1.1  Motivation 

A  problem  of  considerable  and  increasing  importance  within  the 
field  of  composite  materials  is  the  fracture  of  laminated  composite 
structures.  The  anisotropy  of  the  material  greatly  complicates  even  the 
simplist  of  problems.  An  example  of  a  complicated  problem  Is  the  edge 
replica  of  Fig.  (1.1);  it  demonstrates  the  crack  types  present  on  the 
free  edge  of  a  [+302/-302]$  tensile  coupon.  From  Fig.  (1.1)  it  is 
apparent  that  a  crack  can  start  out  in  a  transverse  mode  and  turn  into  a 
delamination  within  Its  growing  length.  Complexities  such  as  this 
require  in-depth  fracture  mechanics  models  which  not  only  predict  at 
what  load  a  crack  will  extend  but  also  the  direction  of  crack  exten¬ 
sion.  This  study  was  undertaken  In  an  attempt  to  develop  a  model  with 
the  capability  of  describing  the  characteristics  of  crack  growth  In 
composites. 

1.2  Literature  Review 

Smith  [1]  has  discussed  limitations  of  some  of  the  current  analyti¬ 
cal  models  for  predicting  crack  growth  characteristics  in  composite 
materials.  Likewise,  most  of  the  previous  finite  element  models  either 
Involve  complex  computational  procedures  or  suffer  from  serious  limita¬ 
tions.  Some  models  distinguish  between  fiber  and  matrix  [2,  3.1,  others 
use  hybrid  or  singular  finite  elements  [4,  5],  and  still  others  assume  a 
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direction  of  crack  growth  [6],  Finite  element  analyses  which  require  a 
distinction  between  the  fiber  and  matrix  materials  are  impractical  for 
most  real  life  situations  since  the  number  of  degrees  of  freedom 
required  in  such  an  analysis  could  easily  exceed  the  capacity  of  most 
computers  currently  available.  While  the  use  of  singular  or  hybrid 
finite  elements  may  give  a  good  representation  of  the  stresses,  strains 
and  displacements  near  a  crack  tip,  they  generally  require  more  computa¬ 
tion  time  and  they  are  not  well  suited  to  explaining  large  scale  crack 
growth  characteristics.  Experimental  results  [7,  8]  have  shown  that  the 
direction  of  crack  extension  can  change  significantly  during  crack 
growth  in  a  laminated  composite  material.  These  results  limit  analyses 
which  assume  a  direction  of  crack  growth  to  specialized  cases  or  to 
small  increments  of  crack  growth. 

There  are  many  different  finite  element  models  that  do  not  require 
specialized  elements.  Two  such  models  are  the  failed  element  approach 
of  Adams  [3],  and  the  modified  crack  closure  approach  of  Rybicki  and 
Kannien  [9],  The  failed  element  approach  assumes  that  when  an  element 
in  an  area  of  high  stress  exhausts  its  strain  energy  capacity,  it 
fails.  From  this,  it  is  assumed  that  a  "crack"  has  formed  and  has  the 
dimensions  of  the  failed  element.  This  approach  has  two  implications, 
the  most  Important  of  which  is  that  a  finite  amount  of  material  is 
removed  from  the  system,  which  in  an  actual  material  is  not  the  case. 

The  other  is  that  the  crack  is  not  likely  to  close  up  on  itself  in 
subsequent  loading  because  of  its  exaggerated  width.  The  modified  crack 
closure  technique  is  based  on  the  crack  closure  integral  and  can  be  used 
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within  the  framework  of  a  linear  elastic  analysis  with  a  relatively 
coarse  mesh.  In  the  modified  crack  closure  technique,  the  crack  closure 
integral  is  evaluated  directly  from  the  nodal  forces  and  displacements 
required  to  close  a  virtual  crack  of  extension,  Aa.  The  modified  crack 
closure  technique  also  has  the  advantage  of  obtaining  mode-I,  mode-II, 
and  mode-I I I  results  in  a  single  analysis. 

Few  fracture  theories  predict  the  direction  of  crack  extension  [2], 
as  well  as  the  external  load  level  which  causes  crack  extension.  Hashin 
[in]  suggested  that  a  failure  criterion  could  be  constructed  which  would 
include  the  plane  on  which  failure  would  occur.  Some  of  the  many  frac- 
ture/fai lure  criteria  which  have  been  used  to  predict  fracture/failure 
of  composite  materials  include  the  Si h  strain  energy  density  criterion 
[Z],  the  Tsai-Wu  failure  criterion  [11],  the  Whitney-Nuismer  point 
stress  criterion  [12],  and  the  Hashin  failure  criterion  [10].  Of  the 
criteria  listed,  only  the  strain  energy  density  criterion  [2]  and  the 
point  stress  criterion  [12]  are  readily  capable  of  predicting  the  direc¬ 
tion  of  crack  extension.  However,  without  modification,  their  use  is 
limited  to  special  cases*  >  I 


1.3  Purpose  of  the  Present  Study 

\ 

The  purpose  of  this  study  was  to  develop  a  finite  element  model 
capable  of  predicting  crack  growth  characteristics  in  composite  materi¬ 
als.  It  was  desired  to  develop  a  model  which  not  only  could  determine 
what  applied  lead  level  would  cause  crack  extension  but  one  which  could 
also  determine  the  direction  of  crack  extension. 


I  > 
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1.4  Basic  Assumptions 

Unless  otherwise  stated,  the  model  developed  was  based  on  the 
following  assumptions: 

(i)  Linear  elastic,  homogeneous  isotropic  or 
homogenous  orthotropic  fibrous  composites 

(ii)  small  displacement  theory 

(iii)  a  crack  extends  from  one  end  only  -  one  crack  at 
a  time 

(iv)  no  variation  of  geometry  in  one  of  the  coordinate 
directions  (i.e.,  plane  stress,  plane  strain  and 
generalized  plane  strain  problems). 

1.5  Description  of  the  Finite  Element  Model 

The  finite  element  model  developed  uses  a  two  dimensional  mesh  of 
four  node,  linear,  isoparametric  elements.  The  model  has  the  capabili¬ 
ties  of  obtaining  either  plane  strain,  plane  stress  or  generalized  plane 
strain  solutions.  (Refer  to  Appendix  A  for  a  description  of  the  finite 
element  method  as  it  applies  to  this  study.)  The  material  models  avail¬ 
able  include  isotropic,  orthotropic  and  laminated  orthotropic,  (off- 
axis),  materials.  (See  Appendix  B  for  an  explanation  of  the  constitu¬ 
tive  relations  for  the  respective  material  models.) 

The  approach  taken  to  the  solution  of  the  crack  problem  was  to 
separate  the  analysis  into  two  main  parts.  In  the  first  part  of  the 
solution  the  direction  of  crack  growth  was  determined  and  in  the  second 
part  the  load  level  which  would  cause  crack  extension  was  determined. 
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The  crack  growth  direction  was  determined  through  the  use  of  sev¬ 
eral  fracture/failure  theories.  The  theories  considered  include  a 
modified  version  of  the  Griffith  criterion  [133.  the  Sih  strain  energy 
density  criterion  [2],  the  Tsai-Wu  failure  criterion  [11],  and  modified 
versions  of  the  Whitney  and  Wuismer  point  stress  theory  [12J,  and  Hashin 
failure  criteria  [10.1.  The  failure  criteria  are  described  in  Chapter  2 
and  the  results  are  compared  with  theory  and  experiment  in  Chapter  3. 

The  determination  of  the  load  level  which  would  cause  crack  exten¬ 
sion  was  made  through  the  use  of  the  modified  crack  closure  method 

[91.  The  modified  crack  closure  method,  as  it  applied  to  this  study,  is 
presented  In  Chapter  2. 

The  main  reasons  behind  the  choice  of  this  solution  approach  was 
that  it  could  he  used  with  a  linear  elastic  analysis,  that  it  could  be 

used  with  a  relatively  coarse  mesh  and  that  it  required  a  minimum  of 
computer  time. 


I*6  Problems  Considered 

The  problems  considered  in  this  study  were: 

(i)  A  mode-I  crack  in  an  infinite  isotropic  plate  with  remote 
loading  of  a.  Fig.  (1.2a),  was  considered  as  a  test  of 
the  crack  closure  technique* 

(ii)  Mixed  mode  cracks  in  infinite  plates  of  isotropic  materials 
Fig.  (1.2b),  were  analyzed  as  a  test  of  the  crack  growth 
direction,  -e^  for  various  angles  of  crack  inclination,  6. 


b.  Mixed  Mode  Crack  In  An  Infinite  Plate 


Figure  1.2  Isotropic  Crack  Problems  Considered 
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(iii)  Off-axis,  unidirectional,  orthotropic  tensile  specimens. 

Fig.  (1.3a),  were  analyzed  for  the  direction  of  crack 
extension  in  fibrous  composites  and  the  results  were  compared 
with  available  experimental  results. 

(i v)  Transverse  cracks  were  introduced  on  the  free  edge  of  a 

laminated  composite  tension  specimen.  Fig.  (1.3b),  and  the 
predicted  crack  paths  were  compared  against  available 
experimental  results. 
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a.  Crack  In  A  Unidirectional 
Tensile  Specimen 


b.  Transverse  Crack  In  A  Laminated 
Tensile  Specimen 


Figure  1.3  Composite  Crack  Problems 
Considered 


Chapter  2 

THEORETICAL  BACKGROUND 

2.1  Criterion  for  Crack  Extension 

2.1.1  Energy  Release  Rate  Concept 

In  the  Griffith  theory  of  classical  fracture  mechanics  [13],  it  is 
assumed  that  strain  energy  is  released  when  a  crack  surface  is  created 
in  a  stressed  body.  The  rate  of  energy  release  when  a  crack  extends 
stably  in  a  body  is  known  as  the  critical  energy  release  rate,  Gc.  The 
critical  energy  release  rate,  Gc,  can  be  determined  experimentally  by  a 
procedure  which  allows  for  stable,  slow  crack  extension,  (see  Ref.  [17] 
for  isotropic  materials  and  [6]  for  composites). 

For  structures,  such  as  a  composite  laminate,  an  existing  crack  may 
or  may  not  grow  under  a  given  state  of  stress.  To  determine  whether  or 
not  an  existing  crack  will  extend,  it  is  necessary  to  calculate  the 
available  energy  release  rate,  G(a),  associated  with  a  crack  of  length 
a.  If  the  available  energy  release  rate,  G(a)j  is  equal  to  the  critical 
energy  release  rate,  Gc  the  crack  will  grow  in  a  stable  fashion.  If  the 
available  energy  release  rate,  G(a),  is  greater  than  the  critical  energy 
release  rate,  Gc,  the  crack  grows  unstably  and  if  the  available  energy 
release  rate,  G(a),  is  less  than  the  critical  energy  release  rate,  Gc, 
the  crack  does  not  extend.  Similarly,  the  external  load  which  first 
causes  the  available  energy  release  rate,  G(a),  to  reach  the  critical 
energy  release  rate  value,  Gc,  is  the  critical  load  and  loads  greater 
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than  the  critical  load  level  result  in  unstable  growth  while  loads  less 
than  the  critical  load  level  do  not  extend  the  crack. 

Mathematically,  the  energy  release  rate,  6(a),  for  a  given  crack  of 
initial  length,  a,  is  defined  as  the  difference  in  the  total  strain 
energy  of  the  structure,  AU,  before  and  after  a  small  crack  extension, 
Aa,  is  introduced,  that  is, 


6(a) 


lim  AU 
Aa+0  Aa 


(2.1) 


2.1.2  Modified  Crack  Closure  Approach 

Irwin  [14]  contended  that  if  a  crack  extends  by  a  small 
amount,  Aa,  the  energy  released  in  the  process  is  equal  to  the  work 
required  to  close  the  crack  back  to  its  original  length.  This  statement 
in  equation  form  is 


G(a) 


1  im  1 
Aa+0  2Aa 


Aa 

J  o«Au  da 
0 


(2.2) 


where  o  Is  the  surface  traction  vector  and  Au  the  displacement  vector 
required  to  close  the  crack  back  to  its  original  length. 

The  modified  crack  closure  technique  of  Rybicki  and  Kannien  [9] 
enables  the  direct  evaluation  of  the  crack  closure  integral  (2.2)  and 
thus  the  energy  release  rate  through  the  use  of  a  finite  element 
model.  The  finite  element  model  starts  with  the  presence  of  an  initial 
crack  of  length,  a.  Fig.  (2.1a),  with  tip  at  node  K.  The  finite  element 
solution  determines  the  displacement  components,  u^,  (where  u^  =>(uk. 
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Vk,  Wk)  of  the  crack  tip  node  K.  An  incremental  crack  extension  Aa  is 
introduced  by  replacing  the  crack  tip  node  K  with  two  separate  nodes  K’ 
and  K11  as  shown  in  Fig.  (2.1b).  With  this  new  crack  geometry  taken 
into  account,  the  finite  element  solution  for  the  nodal  displacements 
J  .  and  Jk”  are  found  for  nodes  K1  and  K' 1 ,  respectively,  under  the 
same  load.  The  crack  extension  is  then  closed  by  applying  equal  and 
opposite  forces  at  nodes  K'  and  K' '  such  that  their  common  displacements 
match  the  displacements  found  earlier  for  node  K,  Fig.  (2.1c).  These 
nodal  forces  can  be  described  by 

p  _v  /c  p  f  )  (2.3) 

^k  '  xk’  yk*  zk' 


The  energy  release  rate  is  then  given  by  [91, 

6(.)  =  tFxk(uk'-uk")*Fjk(«k-vk")*Fzk(»k'-wk'V)]/7M  l*-*> 

By  ..solving  the  forces  and  displacements  into  ;a  “crack  coordinate 
system,"  Fig.  <2. 2).  the  respective  fracture  mode  contributions  to  the 

total  energy  release  rate  can  be  determined,  lhat  is, 

\ 

G j (a)  -  |[F2kcos*-Fjks1n*][cose(»k'-wk”)-s1n»(vk'-vk'')])/2ia  (2.5a) 

A 

G„{a)  -  l[Fykcos*fFzks1ne)[cost(vk'-.k'')«l"fWk,-»k"»!(2‘a  <2-5b> 

5„!<s>  *  i 

\ 

i 
! 


(2.5c) 
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where  Gj,  Gjj  and  Gin  are  the  mode-I,  mode-II  and  moee-III  contribu¬ 
tions  to  the  total  energy  release  rate,  G,  respectively,  i.e., 

G(a)  =  Gj(a)  +  Gjj(a)  +  Gjjj(a)  (2.5d) 

2a. 3  Finite  Element  Considerations  in  Computing  the  Energy 
Release  Rate 

As  previously  mentioned,  the  energy  release  rate  requires  the 
evaluation  of  the  nodal  forces  and  displacements  necessary  to  close  a 
crack  of  extended  length,  a  +  a a,  back  to  its  original  length,  a.  The 
needed  displacements  are  directly  obtained  from  finite  element  solutions 
of  the  initial  and  extended  crack  states.  Fig.  (2.1a)  and  (2.1b), 
respectively.  However,  the  calculation  of  the  required  forces  are  not 
as  obvious.  Rybicki  and  Kanninen  [9]  computed  the  forces  by  placing  a 
very  stiff  "spring"  between  nodes  K'  and  K' ' ,  then  computed  the  force 
components  in  the  "spring."  This  procedure  can  lead  to  unnecessary 
approximation  errors.  An  alternative  approach  will  now  be  presented. 
Consider  three  separate  states.  State  Ho.  1  represents  the  loaded 
initial  state.  Fig.  (2.1a),  where  node  K  displaces  (uk,vk,wk).  The 
finite  element  equations  (Appendix  A)  for  State  No.  1  can  be  written  as, 

=  {Fj}  (2.6a) 

where  [Kj]  is  the  global  stiffness  matrix,  (Sj)  is  the  global  displace¬ 
ment  vector  and  { Fj }  is  the  global  force  vector  for  State  No.  1.  Simi- 
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larly.  State  No.  2,  Fig.  (2.1b),  represents  the  loaded  extended  state 
which  can  be  expressed  by, 

[K2]{62}  «  t^l  (2‘6b) 

and  State  No.  3,  Fig.  (2.1c),  represents  the  loaded  separated  state  with 
applied  forces  to  hold  nodes  K'  and  K1 1  together  -  which  may  be  written 
as, 

[K3]{63}  =  {F3}  (2‘6C) 

Since  the  forces  required  to  hold  nodes  K1  and  !'**  together  are  con- 
tained  within  the  {F3)  vector  it  is  necessary  to  compute  {F3J  .  Since 
the  separated  state  with  applied  forces.  State  3,  Fig.  (2.1c),  is  con¬ 
strained  to  displace  identical  to  that  of  the  unseparated  state.  State 
1,  Fig.  (2.1a),  the  displacement  vector  {6^}  is  the  same  as  {63} 
the  exception  of  the  additional  degrees  of  freedom  (i.e.,  u,  v  and  w 

displacements)  for  the  new  node  created  by  separating  the  crack  tip  node 

v  \ 

into  two  nodes.  Now,  if  {6  '}  is  defined  as  being  the  {Sj}  vector  with 
the  additional  degrees  of  freedom  it  follows  that 

X  \ 

t«3l  *  lv>  \ 

\ 

I 

where  the  additional  degrees  of  freedom  are  specified  as  being  the  same 
as  for  the  initial  crack  tip  node  of  State  1  since  State  3  requires  that 

i 

| 

I 
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the  displacements  of  the  separated  nodes  match  those  of  the  unseparated 
state,  State  1.  Note  that  the  new  node  created  by  separating  the  crack 
tip  node  was  numbered  as  being  the  crack  tip  node  number  plus  one  and 
all  node  numbers  greater  than  the  crack  tip  node  were  re-numbered  as 
being  one  plus  the  node  numbers  that  they  had  in  State  1.  This  re¬ 
numbering  procedure  guarantees  that  the  half  bandwidth  will  „ot  increase 

by  any  more  than  2  for  plane  stress  or  plane  strain  and  by  an,  more  than 
3  for  generalized  plane  strain. 

The  undeformed  mesh  of  the  separated  state.  State  2,  Fig.  (2.lb)  ,s 
identical  to  that  of  State  3,  Fig.  (2.1c),  an<1  stnce  the  stm„ess 

matrices  do  not  change  for  different  loading  conditions,  it  follows  that 


ny  = 


(2.6e) 


Substituting  Eqns.  (2.6d)  and  (2.6e)  into  (2.6c),  the  solution 


for  {F3}  is  found  to  be 


£k23fy}  =  {y  (2.7) 

Hence,  for  a  growing  crack  problem,  the  forces  necessary  to  close  the 
current  crack  extension  are  found  by  simply  multiplying  the  current 
stiffness  matrix  by  the,  modified,  previous  displacement  vector  -  with- 
out  the  addition  of  extra  steps  or  the  introduction  of  unnecessary 
approximation  errors.  Note  that  it  is  not  necessary  to  store  the  entire 
[k2]  stiffness  matrix.  Since  [k2]  is  a  banded  matrix,  the  only  contri- 
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butions  to  the  force  necessary  to  close  the  crack  will  come  from  the 
elements  containing  the  crack  tip  node  K.  Hence,  the  elemental  contri¬ 
butions  to  the  force  at  node  K  can  simply  be  summed  up  to  give  the  total 
force  at  node  K. 

?'2  Criteria  for  Predicting  the  Direction  of  Crack  Extension 
?.2A  Modified  Griffith  Criterion 

The  Griffith  or  energy  release  rate  criterion  states  that  a  crack 
will  extend  when  the  available  energy  release  rate,  G(a),  reaches  or 
exceeds  the  critical  energy  release  rate,  Gc  [13].  In  a  crack  problem 
where  the  crack  extension  direction  is  unknown  the  criterion  should  be 
modified  to  state  that  crack  extension  will  occur  in  the  direction  in 
which  the  available  energy  release  rate,  G(a),  first  reaches  the  criti¬ 
cal  energy  release  rate,  G£. 

If  the  critical  energy  release  rate  is  assumed  independent  of 
direction  then  the  direction  of  crack  extension  can  be  taken  as  the 

j 

direction  of  maximum  available  energy  release  rate  since  this  would  be 

the  direction  which  would  first  reach  or  exceed  the  critical  energy 
release  rate.  \ 

Two  serious  limitations  of  the  modified  Griffith  criterion,  as 
defined  above,  are  that  the  critical  energy  release  rate  may  have  a 
dependence  on  the  mode  of  fracture  in  isotropic  problems  and  it  also 
depends  on  which  direction,  relative  to  the  material  principal  coordi¬ 
nates,  the  crack  extends  in  fibrous  materials.  As  an  example  of  the 

latter  of  the  two  limitations,  consider  two  different  mode-I  cracks 

| 

|  *  ^ 

I 

; 
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extending  in  an  anisotropic  material.  The  first  crack.  Fig.  (2.3a), 
represents  a  mode-I  crack  growing  parallel  to  the  fibers  and  the  second 
crack.  Fig.  (2.3b),  represents  a  mode-I  crack  growing  perpendicular  to 
the  fibers.  Based  on  surface  energy  considerations  [13],  an  approximate 
relation  for  the  mode-I  critical  energy  release  rate,  Gjq,  for  an  iso¬ 
tropic  material  is  given  by 

2 

iro  a 

r  _  **■  /o  o\ 


where  is  the  critical  applied  stress  required  to  cause  crack  exten¬ 
sion  and  E  is  Young's  modulus.  Substituting  the  ultimate  strength  of  a 
composite  material,  (T300/5208  graphite-epoxy),  in  the  transverse  direc¬ 
tion,  Yj,  and  the  modulus  in  the  transverse  direction,  E2,  from  Appendix 
C,  into  Eqn.  (2.8),  gives  an  approximate  value  for  Gjq  for  extension 
parallel  to  the  fibers.  Fig.  (2.3a),  that  is. 


’IC 


tf(Yt)  a 


81.2a 


(2.9a) 


Similarly,  for  the  crack  growing  perpendicular  to  the  fibers.  Fig. 
(2.3b), 


GIC  = 


ir(XT)2a 


7883a 


(2.9b) 


Hence,  the  critical  energy  release  rate  for  the  crack  growing  perpen¬ 
dicular  to  the  fibers  is  roughly  two  orders  of  magnitude  greater  than 


x 
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that  of  the  crack  growing  parallel  to  the  fibers. 

Another  disadvantage  of  the  modified  Griffith  criterion  is  that  it 
requires  an  additional  finite  element  solution  for  each  possible  direc¬ 
tion  that  the  crack  can  extend.  Referring  to  Fig.  (2.4),  for  the  crack 
defined  by  oab  there  are  seven  possible  directions  of  crack  extension, 
from  o  to  c,  from  o  to  d,  from  o  to  e,  from  o  to  f,  from  o  to  g,  from  o 

to  h  and  from  o  to  i.  In  order  to  use  the  modified  Griffith  criterion, 

seven  independent  finite  element  solutions  would  be  required  to  compute 
the  seven  possible  energy  release  rates.  This  is  obviously  time  consum¬ 
ing  and  thus  a  costly  procedure. 

2.2.2  Sih  Strain  Energy  Density  Criterion 

The  strain  energy  density  criterion  [2,  15]  is  based  on  the  local 

value  of  strain  energy  density  in  the  vicinity  of  a  crack  tip,  which  is 

direction  sensitive.  Crack  extension  is  postulated  to  occur  in  the 
direction  of  minimum  strain  energy  density  when  the  strain  energy 
density  factor,  S,  (to  be  defined),  attains  a  critical  value,  Sc. 

For  a  planar  crack  in  an  isotropic  material  under  plane  strain. 
Fig.  (2.5),  the  strain  energy  density  in  the  vicinity  of  the  crack 
tip,  -gy-,  is  given  as,  [15]  \ 

\ 

W  =  ^allkI+2a12kIkII+a22kII+a33kIII^  \ 

+  non-singular  terms 


(2.10) 


u  ORIGINAL  page  is 
OF  POOR  QUALITY 


Finite  Element  Representation 
Of  A  Crack  With  Seven  Possibl 
Extension  Directions 
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where  kj,  kjj  and  kjjj  are  the  mode-I,  mode-II  and  mode-III  stress 
intensity  factors,  respectively. 


all  "  ■y^jC(3-4v-cos<J> ) (l+cos<f> )1 
a12  =  TEG  *  2sin*[cos*-(l-2v)] 
a22  =  igGC4(1-v)(1-c°S(f))+(l+cos<}.)(3cos<t>-l)3 


(2.11a) 

(2.11b) 

(2.11c) 


(2. lid) 


and  G  is  the  shear  modulus  of  elasticity  and  v  is  Possion's  ratio. 

Eqn.  (2.10)  demonstrates  that  the  strain  energy  density  function 
possesses  a  (1/r)  singularity  at  the  crack  tip.  Hence  the  expression 


2  2  2 

S  =  (a  j+2a  I '*'^22^1  I*a33^I  1 1  (2.12) 

i 

I 

! 

represents  the  intensity  of  the  strain  energy  density  field  in  the 
vicinity  of  the  crack  tip.  The  fundamental  hypotheses  on  crack  growth 
in  the  Si h  theory  are  as  follows: 

\ 

(1)  Crack  initiation  takes  place  in  a  direction  determined  by  the 

\ 

stationary  value  of  the  strain-energy  density  factor,  i.e., 

H  *  ».  at  ♦  ■  *0  ;  (2-13) 


/-  '  -- 
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(2)  Crack  extension  occurs  when  the  strain-energy  density  factor 
reaches  a  critical  value,  i.e.. 


•>£  —  S(kj,  kj j,  for  f  —  ^  (2.14) 

Exact  evaluation  of  the  stresses  and  strains  in  the  vicinity  of  the 
crack  tip  with  the  current  finite  element  model  is  not  guaranteed  since 
there  exists  a  geometric  singularity  at  the  crack  tip  which  cannot  be 
accurately  modeled  with  the  linear-elastic  analysis  formulated  herein. 
Hence,  it  is  not  possible  to  use  the  Sih  theory  to  determine  when  the 
crack  will  extend.  However,  the  Sih  theory  can  be  used  to  determine  the 
direction  of  crack  propagation  in  isotropic  materials. 

From  continuum  mechanics  [16]  it  is  possible  to  write  an  alterna¬ 
tive  form  of  the  strain  energy  density  at  a  point  in  a  stressed  body, 
i.e.. 


dU  1. 

dV  "  2^axxexx+0yyeyy+O2zezz+Tyz1ryz+TxzYxz+TjcyTxy^ 


(2.15) 


Neglecting  the  non-singular  terms  in  Eqn,  (2.10)  and  substituting  In  the 
expression  for  the  strain  energy  density  factor  of  Eqn.  (2.12)  gives. 


dU  _  S 
q'V  ~  r 


(2.16) 


/ 
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Combining  Eqns.  (2.15)  and  (2.16),  an  alternative  equation  for  the 
strain  energy  density  factor,  S,  is  found  to  be: 


£(a 


xxexx+ayyeyy+  azzezz+TyzYyz+TxzYxz+TxyYxy' 


(2.17) 


There  is  one  serious  limitation  to  the  use  of  the  strain  energy 
density  theory  in  the  current  study.  This  limitation  is  that  the  theory 
does  not  account  for  the  anisotropic  strength  characteristics  of  the 
material.  Since  such  properties  must  be  accounted  for  in  fibrous  com¬ 
posites,  the  Si h  strain  energy  density  theory  is  limited  to  isotropic 
material  applications.  It  should  be  noted  that  the  strain  energy 
density  theory  has  been  used  in  the  past  to  predict  crack  growth  charac¬ 
teristics  in  composite  materials  [2].  However,  the  success  of  such 
studies  resulted  from  assuming  that  the  crack  was  situated  entirely 
within  the  Isotropic  matrix  between  fibers. 

The  procedure  for  implementing  the  strain  energy  density  criterion 
in  the  finite  element  model  is  briefly  described  as  follows:  First,  the 
possible  directions  of  crack  extension  in  the  model  are  identified  by 
the  element  sides  containing  the  crack  tip  node,  (node-0  in  Fig. 

(2.4)).  Second,  the  stresses  and  strains  are  calculated  in  the  adjoin¬ 
ing  elements  at  the  element  corners,  (points  c  thru  i  in  Fig.  (2.4)). 
Third,  Eqn.  (2.17)  is  used  to  calculate  the  strain  energy  density  fac¬ 
tor,  S,  at  the  respective  points.  Last,  the  crack  is  assumed  to  grow  in 
the  direction  in  which  S  is  a  minimum. 
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2.2.3  Tsai-Wu  Failure  Criterion 

Tsai  and  Wu  [11]  postulated  that  a  failure  surface  in  stress  space 
exists  in  the  form: 


F.s  .+F .  .s  .s  .  =  1 
ij  i  J 


i,  j  =  I. 


(2.18) 


where  Fi  and  Fi  .•  are  strength  tensors  of  second  and  fourth  order, 
respectively,  and  s^  represent  a  contracted  form  of  the  stress  tensor 
components  in  material  principal  coordinates.  For  an  orthotropic  lamina 
under  plane  stress  conditions,  Eqn.  (2.18)  becomes: 


(xI+Xi)sll+  (y;  +  T1')s22“(X'X  )S11"(Y  Y  )s22+(  2  12 

T  c  T  c  I  c  i  c 


+2F12SUS22  =  1 


(2.19) 


where  Xy  and  Yy  represent  the  tensile  strength  of  the  material  in  the 
fiber  and  transverse  directions,  respectively,  Xc  and  Yc  represent  the 
compressive  strengths,  Sj2  represents  the  shear  strength  in  the  1-2 
plane  and  F12  is  an  interaction  term  which  must  be  determined  from  a 

biaxial  strength  test.  \ 

As  a  failure  theory,  the  Tsai-Wu  criterion  has  several  advanta¬ 
ges.  These  advantages  include,  (1)  invariance  under  rotation  of  coordi 
nates,  (2)  transformation  via  known  tensor  transformation  laws;  and  (3) 
symmetry  properties  akin  to  those  of  the  stiffnesses j and  compliances. 
However,  for  use  in  this  study,  it  has  two  serious  limitations.  The 


J 


28 

first  is  that  it  fails  to  account  for  differences  in  creating  new  frac¬ 
ture  surfaces  at  various  angles  to  the  material  principal  coordinates 
and  second  is  that  it  makes  no  distinction  between  tensile  and  compres¬ 
sive  failure  which  could  cause  the  Tsai-Wu  criterion  to  choose  a  com¬ 
pressive  direction  of  crack  propagation. 

The  Tsai-Wu  criterion  was  incorporated  into  the  finite  element 
model  by  assuming  that  the  crack  would  extend  in  the  direction  where  the 
value  of  the  Tsai-Wu  polynomial  reached  a  maximum.  Referring  to  Fig. 

(2.4) ,  the  stresses,  a..,  are  calculated  at  some  fixed  distance,  rQ, 
away  from  the  crack  tip  at  the  various  locations  dictated  by  the  element 
sides  incorporating  the  crack  tip  node,  (points  c  thru  i  in  Fig. 

(2.4) ).  Next,  the  values  cf  the  Tsai-Wu  polynomial  were  computed  at 
these  points  through  the  use  of  Eqn.  (2.19),  (for  the  case  of  orthotro¬ 
pic  plane  stress).  Last,  the  crack  was  assumed  to  extend  in  the  direc¬ 
tion  for  which  the  Tsai-Wu  polynomial  reached  a  maximum. 

The  choice  of  r0  is  arbitrary  within  certain  limitations.  These 
limitations  are  that  rQ  should  be  greater  than  zero  and  less  than  the 
longest  possible  path  that  the  crack  extension  could  take  and  not  extend 
through  more  than  one  element.  Referring  to  Fig.  (2.4),  if  the  possible 
path  of  crack  extension  were  as  shown  from  node  o  to  c,  o  to  d,  o  to  e, 
o  to  f ,  o  to  g,  o  to  h  or  o  to  i ,  then  would  be  limited  to  that  of 
the  segment  from  o  to  g  since  that  is  the  direction  of  longest  possible 
single  element  crack  extension.  Note  that  it  would  be  impractical  to 
use  a  distance  greater  than  og  since  this  would  require  using  stresses 
from  an  element  outside  those  adjacent  to  the  crack  tip. 
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2.2.4  Modified  Point  Stress  and  Hashin  Criteria 

The  criteria  for  predicting  the  crack  growth  direction  considered 
up  to  this  point,  in  their  present  form,  are  all  unsatisfactory  for 
anisotropic  materials.  As  pointed  out  in  Sections  2. 2. 1-3,  they  all 
fail  to  account  for  differences  in  the  work  required  to  create  a  new 
area  of  crack  surface  at  different  directions  In  an  anisotropic  mate¬ 
rial.  Since  it  is  imperative  that  this  distinction  be  accounted  for  in 
this  study,  two  new  criteria  are  proposed.  The  first  criterion  consid¬ 
ered  is  a  modification  of  the  point  stress  criterion  of  Whitney  and 
Nuismer  [12],  and  fcne  second  is  a  modification  of  the  Hashin  criterion 
[10]. 

2. 2. 4.1  Modified  Point  Stress  Criterion 

The  point  stress  criterion  of  Whitney  and  Nuismer  [12]  assumes  that 
failure  of  a  notched  laminate  occurs  when  the  local  stress  at  a  certain 
distance,  r  ,  from  the  notch  tip  reaches  the  strength  of  toe  unnotched 
laminate. 

The  modified  point  stress  criterion  of  this  study  assumes  that  a 

crack  will  grow  in  the  direction  of  the  maximum  ratio  of  normal  stress 

to  strength  at  a  certain  distance,  r  ,  from  the  tip  of  an  existing 

crack.  Note  that  this  is  equivalent  to  assuming  that  a  crack  will  grow 

perpendicular  to  the  plane  of  maximum  tensile  stress  in  an  isotropic 

material.  Referring  to  Fig.  (2. A),  the  normal  stress,  (r,<|>),  is 

99 

calculated  at  some  fixed  distance,  r  ,  away  from  the  crack  tip.  The 
stress,  o  (r  ,<(>),  is  then  divided  by  the  tensile  strength  of  the  mate- 

99  0 
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Figure  2.6  Modified  Point  Stress  Parameters 
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rial,  Tx  (to  be  defined),  normal  to  the  direction  of  crack  exten- 

99 

sion.  This  ratio,  R(rQ,9),  is  then  used  to  predict  the  direction  of 
crack  extension  by  assuming  that  the  crack  will  extend  in  the  direction 
for  which  R(rQ,9)  reaches  a  maximum.  Quantitatively,  the  ratio  is 
defined  as: 

R(r  ,♦)  .  °h(,'0’',)  (2.20) 

V*> 

The  value  of  rQ  used  in  this  criterion  is  subject  to  the  same 

limitations  as  the  value  of  r  in  the  Tsai-Wu  criteria.  Section  2.2.3. 

o 

That  is,  rQ  should  he  greater  than  zero  and  less  than  or  equal  to  the 
longest  path  of  possible  crack  extension  while  not  extending  through 
more  than  one  finite  element. 

2. 2. 4. 2  Strength,  T  ,  Along  a  Given  Plane  in  Anisotropic  Materials 
99 

The  strength,  T  (9),  normal  to  a  given  direction,  is  taken  as  the 
99 

normal  stress  required  to  fail  an  infinitesimal  element  of  anisotropic 
material  along  a  qiven  plane.  In  the  finite  element  solution  the  finite 
element  sides  dictate  the  directions  of  possible  crack  extension. 

T..  was  taken  in  this  form  to  account  for  differences  in  the  energy 

99 

required  to  create  new  crack  surfaces  at  arbitrary  angles  with  the 
material  principal  coordinate  system.  Such  a  definition  is  necessary  in 
the  proposed  model  to  permit  selection  of  the  proper  direction  of  crack 
growth.  Further,  as  shown  by  Herakovich  [7],  for  example,  failure  of 
individual  lamina  in  a  laminate  can  occur  along  planes  which  are  neither 
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parallel  to  nor  perpendicular  to  the  fibers.  It  appears  to  be  most 

difficult  If  not  Impossible  to  experimentally  verify  a  relation 

for  I  U)  by  testing  unidirectional  laminates  since  only  one  direction 
♦  ♦ 

of  crack  extension  would  be  present  for  a  given  laminate.  However, 
three  obvious  conditions  must  he  met  by  such  an  expression.  The  condi¬ 
tions  are: 

\.  Tor  an  Isotropic  material,  the  strength.  T^(*).  should  he  constant 
and  equal  to  the  ultimate  strength  of  the  material,  °u»  Independent 

of  4>.  ' 

?.  For  a  crack  extending  parallel  to  the  fibers  In  a  unidirectional 
composite.  T  U)  should  be  equal  to  the  transverse  tensile 
strength  of  the  material.  Yy. 

3.  For  a  crack  extending  perpendicular  to  the  fibers*  ^^(4*)  should  be 
equal  to  the  tensile  strength  of  the  material  In  the  fiber  dlrec- 

tlon,  Ky. 

For  a  unidirectional  laminate  under  plane  stress,  as  shown  In  Fig. 
(?.7).  a  simple  relationship  for  T^(*)  can  t?e  postulated.  Removing  an 
infinitesimal  element  at  (ro.»).  Fig.  (?.R«K  and  defining  the  angle  d 
as  the  difference  between  the  fiber  angle.  0.  and  the  assumed  crack 
extension  angle.  0,  that  Is 

\ 

'  ft  e  0  -  4  \  (?•?!) 

i 

i 

It  Is  then  possible  to  Isolate  yet  another  Infinitesimal  element.  Fig. 
(?.ah).  which  gives  the  orientation  of  the  crack  extension  relative  to 
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the  material  principal  directions.  Then,  assuming  that  the  strength, 

»  ’s  balanced  by  only  the  transverse  and  longitudinal  strength  of 
the  material  Yy  and  Xy,  respectively,  the  relation  for  T  ($)  is 
obtained  by  summing  forces  in  the  4>  direction  as. 


th(«)  =  xrsin  0  + 


YyCOS26 


(plane  stress) 


(2.22) 


Testing  Eqn.  (2.22)  against  the  three  conditions. 

Isotropic  material:  X^.  =  Y^  = 

2  2 

T..  =  o  sin  6  +  o  cos  B  =  a 

99  U  U  U 

Fracture  parallel  to  fibers:  B  =  0° 

T  =  Xysin2(0°)  +  Yycos2(0°)  =  Yy 

Fracture  perpendicular  to  fibers:  B  =  90° 

Tw  =  Xysin2(90°)  +  YTcos2(90°)  =  Xy 

I 

The  three  conditions  specified  herein  are  satisfied.  Hence,  in  princi¬ 
ple,  the  expression  is  acceptable.  A  plot  of  Txi  vs.  <J>  for  plane  stress 

9<j> 

is  shown  in  Fig.  (2.9)  for  various  values  of  0  for  T300/5208  graphite- 
epoxy,  (properties  from  Appendix  C).  The  maximumsN in  Fig.  (2.9)  repre¬ 
sent  the  combinations  of  angles  for  which  a  crack  grows  perpendicular  to 
fibers  and  the  minimums  represent  crack  growth  parallel  to  fibers. 


(2.23a) 


(2.23b) 


(2.23c) 


Wksi) 
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For  a  transverse  crack  extending  through  a  laminate,  Fig,  (2.10),  a 
similar  relationship  can  be  derived.  Removing  the  lamina  at  the  crack 
tip  and  rotating  it.  about  the  z'  axis  into  material  principal  coordi¬ 
nates,  Fig.  (2.11a),  gives  the  geometry  necessary  to  compute  the 
strengths  Txx(8),  Tyy  (0 )  and  Tzz (6 )  in  the  x',  y',  and  z*  directions, 
respectively..  Removing  an  infinitesimal  element,  oab,  from  Fig. 
(2,11a),  applying  the  normal  strengths  xT,  YT  and  Txx’  in  the  1',  2'  and 
x'  directions,  respectively.  Fig.  (2. lib),  and  summing  forces  in  the  x'- 

l 

direction  gives  the  relation  for  Txx(0)  as, 

Txx(0)  =  YTsin2e’  +  xjcos20 '  (2.24a) 

Similarly,  removing  element,  ocd,  from  Fig.  (2.11a),  applying 
*  1  « 

strengths  XT,  Yy,  and  Tyy,  Fig.  (2.11c),  then  summing  forces  in  the  y'- 
direction  yields 

Tyy(e )  =  xjsin2e  +  YjCOS20  (2.24b) 

Since  the  z'  and  3’  axis  in  the  lamina  coordinate  system  are  the  same  as 
in  the  laminate  coordinate  system, 

T2Z(6)  =  4  (2.24c) 

I 

The  strength,  T^($),  normal  to  a  free  edge  crack  extension.  Fig. 

(2.12a),  is  found  by  removing  an  infinitesimal  element.  Fig.  (2.12b), 


applying  the  respective  strengths  and  summing  forces  in  the  *  direction 
which  gives 


V*)  »  Vin2*  +  T„cos24 


(2.25) 


Substituting  Eqns.  (2.24a),  (2.24b)  and  (2.24c)  Into  Eqn.  (2.25)  yields 
that 


i  ?  i  i  ?  '  2 


T  (4.)  =  (X-j-si n  9  +YTcos  9  )sin  <t»  +  ZTcos  * 


(2.26) 


Testing  Eqn.  (2.26)  against  the  three  conditions. 


I  t  I 


Isotropic  material:  Xy  =  Yy  =  Zy  =  °u 

I 

'  4"t> 


t'  =  (ous1n29+3ucos2e)sin24>  +  oucos2<j>  =  ou 


(2.27a) 


l 

Fracture  parallel  to  fibers:  9  =0°,  4>  =  90° 

T*  =  (x'sin2(0°)+Y'cos2(0°))sin2(90°)  +  zjcos2(90°)  =  vj  (2.27b) 
66  T  »  l 


Fracture  perpendicular  to  fibers:  9  =90°,  ♦  "  90°  _ 

j'  =  (xlsi n2(90° )+y1cos2(90° )  )sin2(90° )  +  ZyCOS2(90°)  =  Xy 
66  '  T  * 


(2.27c) 


\ 


Again,  the  three  conditions  are  satisfied  so,  in  principle,  the  relation 
is  acceptable.  A  plot  of  vs.  *  for  free  edge; crack  extension  is 
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shown  in  Fig.  (2.13)  for  various  values  of  fiber  orientation,  e  ,  for 
T3D0/5208  graphite-epoxy,  (properties  from  Appendix  C).  As  in  the  plane 

I  I 

stress  case,  the  values  of  $,  (for  9  =  90°),  for  which  is  a  maximum 

represents  crack  extension  perpendicular  to  fibers  while  a  minimum 

I 

T^,  represents  crack  extension  parallel  to  fibers.  The  constant 

t'  for  0*  =  0°  represents  matrix  mode  failure  independent  of  <j>.  Last, 

the  curve  for  e'  =  45°  has  minimums  at  T, .  =  YT,  which  represents  matrix 

99  I 

I 

mode  failure  and  the  maximums  never  reach  Xy  since  some  combination  of 
fibers  and  matrix,  thru  the  width,  is  always  involved  for  the  free  edge 
cracks  considered.  Fig.  (2.10). 


2. 2.4. 3  Modified  Hashin  Criterion 

The  Hashin  failure  criterion  [10]  assumes  that  failure  of  a  trans¬ 
versely  isotropic  material  will  occur  in  a  tensile  fiber 
mode,  Cjj  >  0,  when: 


(— )2  +  —  (°122+0l|)  "  1 


(2.28) 


>12 


where  Xy  is  the  tensile  failure  stress  in  the  fiber  direction  and  Sj,,  is 
the  axial  failure  shear  stress.  The  Hashin  criterion  also  assumes 
failure  to  occur  in  a  tensile  matrix  mode,  +  <?33  >  0,  when: 


2  (°22+<J33^  +  .2  ^°23'°22a33 

YT  b23 


)  +  — (ai2+a13^  s  1 


(2.29) 


12 


(ksi) 
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Figure  2.13  Free  Edge  T^vs.0For  Various  Values  of  Of  0 
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where  and  are  the  transverse  tensile  and  shear  strength  of  the 
material,  respectively.  For  the  case  of  plane  stress  in  the  1-2  plane, 
the  criterion  for  fiber  mode  tensile  failure  becomes: 


°il  ^  °l?~ 

(r1)  +  (-r1)  =  1 

T  12 


and  for  matrix  mode  tensile  failure: 


°22  ^  °12  ^ 

(#)  +(r)  =i 

rT  i2 


(2.30) 


(2.31) 


where  is  the  shear  strength  of  the  material  in  the  1-2  plane. 


The  Hashin  criterion  does  account  for  distinct  differences  in  fiber  and 
matrix  mode  failure.  However,  it  does  not  account  for  an  arbitrary 
combination  of  matrix  and  fiber  mode  failure.  (In  the  actual  failure  of 
composite  laminates  this  feature  is  necessary  since  as  Herakovich  [7], 
for  example,  has  shown  failure  can  occur  on  a  plane  which  is  neither 

parallel  to  nor  perpendicular  to  the  fibers.) 

\  \ 

Hashin  [10]  proposed  that  a  similar  criterion  for  the  failure  of 
composite  materials  could  be  developed  to  include  the  plane  on  which 

failure  occurred,  $  ,  Fig.  (2.14).  Such  a  criterion  would  predict 

°  \ 
failure  when  some  function  g(o..,  $)  satisfied  the  condition: 

^  J 

<f>)  =  i 


(2.32) 
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Figure  2.14  Composite  Laminate  With  Modified 
Hashin  Parameters 


and  that  failure  would  occur  on  a  plane  defined  by,  for  which  (2.32) 
was  first  satisfied  under  monotonical ly  increasing  load. 

The  modified  Hashin  criterion  of  this  study  assumes  that  a  given 
crack  in  a  composite  material  will  extend  in  the  direction,  4>q,  in  which 
the  left  side  of  Eqn.  (2.32)  reaches  a  maximum  when  evaluated  at  some 
fixed  distance,  rQ,  from  the  crack  tip.  (Note:  r^  is  subject  to  the 
same  limitations  as  in  Sections  2.2.3  and  2. 2. 4.1). 

The  development  of  the  modified  Hashin  criterion  as  used  in  this 
study  is  based  on  developing  an  expression  for  Eqn.  (2.32).  As  in  the 
case  of  the  point  stress  criterion,  it  appears  impossible  to  test  such 
an  expression  experimentally.  However,  two  obvious  conditions  should  be 
met  by  such  a  criterion,  they  are: 

1)  For  a  crack  extending  parallel  to  the  fibers,  the  criteria  should 
give  back  the  Hashin  criterion  for  tensile  matrix  mode  failure, 

Eqn.  (2.31). 

I 

2)  For  a  crack  extending  perpendicular  to'  the  fibers,  the  criteria 
should  give  back  the  Hashin  criterion  for  tensile  fiber  mode  fail¬ 
ure,  Eqn.  (2.30). 

Proceeding  along  the  same  line  as  Hashin  [10],  if  the  failure  criterion 
is  taken  as,  . . \ 

{T^)  +  (r*-)  =  >!  0  (2.33) 

4<J>  nj) 

where  T4,$  and  Tr(},  dre  the  normal  and  shear  strengths  and  and  o  are 
the  normal  and  shear  stresses,  respectively,  on  the  plane  of  crack 


47 


0RI6INAL 

O?  POOH 

extension.  Fig.  (2.14),  then  the  direction  of  crack  extension  is  given 
by  the  value  of  4>  for  which  the  left  side  of  Egn.  (2.33)  reaches  a 
maximum,  provided  that  is  greater  than  zero.  (A  less  than 

zero  would  represent  crack  closure.)  For  the  case  of  plane  stress  of  a 
unidirectional  composite  laminate.  Fig.  (2.7),  the  normal  strength, 

T  (<f>),  was  derived  in  Section  2. 2. 4. 2  and  is  given  by  Fqn.  (2.22).  if 
the  shear  strength,  Tr^(<i>)*  is  simply  taken  as  being  is 

T  =  S10  (plane  stress)  (2.34) 

r<fr  12 

Then  the  failure  criterion  is  complete.  Testing  the  criteria  of  Eqn. 
(2.33)  against  the  two  conditions  specified  herein, 

(i)  Fracture  perpendicular  to  fibers:  4>  =  90  +  6,  equilibrium  of  an 
element.  Fig.  (2.15a)  gives 

2  2  . 
o  =  a  sin  4>  +  a  cos  4>  -  2t  sin$cos4> 
y y  zz 

'  2.2 
°r<j,  =  cos4>si n<f>  (°zz_0yy)  +  Tyz(cos  n  <t>) 

substituting  $  =  90  +  6, 


(2.35a) 

(2.35b) 


PAGE  US 
QUALITY 


oyysin?(90+e)+azzcos2(90+e)-2T^zsin(9049)cos(90+e) 
2  2 

o  cos  sin  6+2t  sin0cos9 

yy  zz  yz 


(2.35c) 
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(ii)  Fracture  parallel  to  fibers:  <t>  =  0  Substituting,  <J»  =  0,  into 
Eqns.  (2.35a)  and  (2.35b)  gives  that: 

2  2 

a..  =  o  sin  6+o_, cos  9-2x  ,sin9cos9 
4|9  yy  zz  yz 

and 

or4i  =  cos6sin9(ozz  cyy)  +  xyz(cos29-sin29) 

but  from  Fig.  (2.15c): 

2  2 

°22  =  °yys1"n  6  +  °zzcos  0  ”  2Tyzsin9cos0 
and 

2  2 

Oj2  =  cos9si n9 (azz-ayy)  +  Tyz(cos  9 -sin  9) 

Comparing  Eqns.  (2.36c)  and  (2.36d)  with  (2.36a)  and  (2.36b),  respec¬ 
tively. 


(2.36c) 

(2.36d) 


(2.36a) 

(2.36b) 


quautv 


and 


(2.36e) 


(2.36f) 


Also,  for  fracture  parallel  to  the  fibers,  Eqn.  (2.22)  reduces  to  that 
of  Eqn.  (2.23b).  Substituting  Eqns.  (2.23b),  (2.34),  (2.36e)  and 
(2.36f)  into  Eqn.  (2.33)  yields. 


°?2  2  °12  ^ 

( y22)  +  (#) 

T  *12 


1 


(2.36g) 
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which  is  precisely  the  Hashin  criterion  for  tensile  matrix  mode  failure, 
Eqn.  (2.31).  Since  the  relation  developed  gives  back  the  Hashin  cri¬ 
teria  for  tensile  fiber  mode  and  tensile  matrix  mode  failure  it  is  in 
principle  acceptable. 

2. 2. 4. 4  Finite  Element  Considerations  in  Implementing  the  Modified 
Point  Stress  and  Kashin  Criteria 

The  steps  in  the  implementation  of  the  modified  point  stress  and 
Hashin  criteria  are  similar  to  those  in  implementing  the  Si h  and  Tsai-Wu 
criteria.  Briefly: 

(i)  Determine  the  elements  containing  the  crack  tip  node.  This  gives 

the  possible  directions  of  crack  extension. 

(ii)  Find  the  minimum  element  side  length.  This  gives  r  . 

o 

(iii)  Compute  the  stresses  along  the  element  sides,  which  define  the 

possible  directions  of  crack  extension,  at  r  . 

o 

(i v )  Use  Eqn.  (2.20)  for  the  point  stress  criterion  or  Eqn.  (2.33)  for 
the  Hashin  criterion  and  assume  crack  extension  in  the  direction 
which  makes  (2.20),  for  point  stress,  or  (2.33),  for  the  Hashin 
criterion,  a  maximum. 


Chapter  3 
RESULTS 


3.1  Isotropic  Cases 

3.1.1  Mode-I  Crack  in  an  Infinite  Plate 

The  classical  problem  of  a  mode-I  crack  in  an  infinite  plate.  Fig. 
(1.2a),  was  run  as  a  test  of  the  energy  release  rate  formulation. 
Results  for  two  different  finite  element  meshes,  one  being  much  finer 
than  the  other,  were  generated  for  comparison.  The  computed  energy 
release  rates  were  converted  to  stress  intensity  factors  for  ease  of 
comparison  with  theory. 

The  boundary  condition  for  the  fine  mesh,  a  306  element  x  338  node 
mesh.  Fig.  (0.1)  considered  a  full  crack  model  assumed  specified  dis¬ 
placement  loading.  Referring  to  Fig.  (3.1a),  the  boundary  conditions 
were: 


at  y  =  0:  v(y  =  0,z)  =  -6 
at  y  =  L:  v(y  =  L,z)  =  6 
at  z  =  0,2B:  traction  free 


(3.1a) 
(3 . lb) 
(3.1c) 


The  input  parameters;  a,  6,  L,  B,  A  and  P  were  taken  as, 

"\ 

a  =  0.5'',  Aa  =  0.2a,  L  =  40a,  B  =  10a  I 

A  =  L/2,  6  =  .002a,  5  =  90°  1  (3.2) 
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Figure  3.1  Finite  Element  Crack  Models 
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The  problem  was  treated  as  a  case  of  plane  strain,  constitutive  equation 
given  by  Eqn.  (B.18),  with  material  properties  listed  in  Appendix  C. 

The  crack  was  assumed  to  be  a  virtual  crack,  i.e.,  having  a  width  equal 
to  zero. 

The  coarse  mesh,  a  68  element  x  82  node  mesh.  Fig.  (D.2),  assumed 
symmetry  about  the  midplane  and  also  used  specified  displacement  load¬ 
ing.  Referring  to  Fig.  (3.1b),  the  boundary  conditions  were: 


at  y  =  0: 

v(y  =  0,z)  =  -6 

(3.3a) 

at  y  =  L: 

v(y  =  L,z)  =  6 

(3.3b) 

at  z  =  0: 

Traction  free 

(3.3d) 

at  z  =  B: 

w(y,  z  =  B)  =  0 

(3.3d) 

The  Input  parameters  were  taken  the  same  as  for  the  fine  mesh,  Eqn. 
(3.2). 

The  theoretical  stress  intensity  factor,  Kj,  for  a  mode-I  crack  in 
an  infinite  plate.  Fig.  (3.1),  is  given  by  [17]  as, 

Kj  =  o(ita)1/2  (3.4) 

and  the  relation  between  the  mode-I  stress  intensity  factor,  Kj,  and  the 
mode-I  energy  release  rate,  Gj ,  for  plane  strain  is  given  by  [17]  as, 

Gf£  i/p 

Ki =  ( — i)  *  (3.5) 

1-v 
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where  E  is  Young's  modulus  and  v  is  Possion's  ratio.  Since  the  theore¬ 
tical  stress  intensity  factor,  Eqn.  (3.4),  is  computed  using  the  remote 
stress,  o,  the  applied  displacement,  6,  must  be  converted  to  an  equiva¬ 
lent  stress.  From  Eqn.  (B.18)  this  relation  is  given  by 


where 


(3.6) 


(3.7) 


The  finite  element  results  for  the  two  meshes  considered  are  pre¬ 
sented  in  Table  (3.1).  The  stresses,  o,  were  computed  using  Eqns.  (3.6) 
and  (3.7).  The  comparisons  of  the  predicted  stress  intensity  factors 
with  the  theoretical  stress  intensity  factors  is  shown  in  Table  (3.2). 
The  results  indicated  pure  mode-I  crack  extension,  i.e.,  -6  =  0°  in 

Fig.  (3.1a),  as  expected.  Note  that  theoretical  values  for  two  differ¬ 
ent  initial  crack  lengths,  a,  are  shown.  The  first  crack  length  repre¬ 
sents  the  actual,  initial,  crack  length  and  the  second  represents  the 
actual,  initial,  crack  length  plus  the  crack  extension.  The  two  crack 
lengths  were  considered  to  demonstrate  the  error  incurred  in  using  a 
finite  crack  extension.  That  is,  the  theoretical  energy  release  rate, 
Eqn.  (2.1),  is  based  on  an  infinitesimal  crack  extension,  i.e.,  limit 
4a  ♦  0,  whereas  the  finite  element  model  introduces  finite  crack 

extension,  4a.  The  results  of  Table  (3.2)  indicated  that  the  fine  mesh 

\ 

gave  better  results  than  the  coarse  mesh  and  that  using  the  crack  length 


Table  3.1  Finite  Element  Energy  Release  Rate  Predictions 
For  A  Mode-1  Crack  In  An  Infinite  Plate 
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plus  the  extension  in  the  theoretical  stress  intensity  factor  compared 
better  with  the  finite  element  results. 

3.1.2  Mixed  Mode  Fracture  of  an  Infinite  Isotropic  Plate 

The  problem  of  a  mixed  mode  crack  in  an  infinite  isotropic  plate. 
Fig.  (1.2b),  was  run  as  a  test  of  the  crack  growth  direction.  Results 
were  compared  with  the  results  predicted  by  the  Sih  strain  energy 
density  theory  [2],  The  mesh  used  consisted  of  306  elements  and  338 
nodes.  Fig.  (E.l). 

The  boundary  conditions  used  were  identical  to  those  of  the  fine 
mesh  for  the  mode-I  crack  in  an  infinite  plate,  i.e.,  Eqns.  (3.1a), 
(3.1b)  and  (3.1c).  Referring  to  Fig.  (3.1a),  the  input  parameters; 
a,  6,  L,  8,  and  A  were  taken  as, 

a=“Mr»  L  =  40a»  B  =  10a 

6=0.001",  A  =  L/2  (3.8) 

i 

t 

while  the  crack  inclination  angle,  &,  was  varied  from  30  to  90  degrees. 

While  simple  relations  for  the  mode-I  stress  intensity  factor,  Kj, 
and  the  mode-II  stress  intensity  factor,  Kjj,  exist  for  the  mixed  mode 
problem  of  Fig.  (1.2b),  reference  [17]  points  out  that  there  is  no  known 

relation  between  the  energy  release  rates  and  stress  intensity  factors 

\ 

for  such  a  problem.  Hence,  discussion  of  the  results  for  the  mixed-mode 
crack  problem  of  Fig.  (1.2b)  will  be  limited  to  the  crack  extension 
direction. 
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The  theoretical  crack  extension  direction,  9q,  was  computed  from 
the  Sih  strain  energy  density  theory  [18],  by  solving  for  6^  from  Eqn. 
(3.9),  below: 

2(l-v)sin(eo-29)  -  {2sin[2(0o-9)]}  -  sin26Q  =  0  (3.9) 

Finite  element  results  were  generated  for  various  crack  inclination 
angles,  B,  through  the  use  of  the  modified  Griffith  criterion.  Section 
2.2.1,  the  Sih  strain  energy  density  theory,  Section  2.2.2,  and  the 
modified  point  stress  and  Hashin  criterion.  Section  2.2.4.  A  plot  of 
the  theoretical  and  finite  element  crack  extension  direction,  6  ,  as  a 
function  of  the  crack  inclination  angle,  B,  is  shown  in  Fig.  (3.2). 

The  results  for  the  modified  Griffith  criterion  were  generated  by 
evaluating  the  value  of  the  crack  closure  integral,  through  the  use  of 
Eqn.  (2.4),  for  all  possible  paths  of  crack  extension  present,  e.g.,  see 
Fig.  (2.4),  and  assuming  that  crack  extension  would  occur  .in  the  direc¬ 
tion  of  a  maximum  energy  release  rate.  The  results  consistently  predic¬ 
ted  crack  extension  in  a  direction  in  which  the  mode-I  energy  release 
rate,  Gj(a),  made  up  99%  or  more  of  the  total  energy  release  rate. 

The  results  for  the  Sih  strain  energy  density  criterion  were  evalu¬ 
ated  using  Eqns.  (2.13)  and  (2.17).  The  strain  energy  density  factor 
was  considered  as  the  sum  of  two  components  [18],  one  due  to  a  change  in 
volume,  Sv,  and  one  due  to  a  change  in  shape,  Sd, 


(3.10) 
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Figure  3.2  Crack  Extension  Angle,— Qq,  vs.  Crack 
Inclination  Angle,  R  ,  For  Isotropic  Mixed  Mode  Fracture 


The  results  always  predicted  crack  extension  in  a  direction  in 

which  S  >  >  S..  This  conforms  with  the  concept  [18]  that  fracture 
v  d 

occurs  along  a  plane  where  Sy  >  S^. 

The  modified  point  stress  cases  were  evaluated  using  Eqns.  (2.20) 

and  (2.22)  with 

*T  *  *T  *  %  <3-U) 

The  crack  was  assumed  to  extend  in  the  direction  for  which,  R(<t> ) ,  Eqn. 
(2.20)  was  a  maximum. 

The  results  for  the  modified  Hashin  criterion  were  generated 
through  the  use  of  Eqn.  (2.33),  with 


(3.12) 


Crack  extension  was  assumed  to  occur  in  the  direction  which  maximized 
the  left  hand  side  of  Eqn.  (2.33). 

All  of  the  tested  theories  gave  identical  results,  as  indicated  in 
Fig.  (3.2).  The  small  differences  between  the  theoretical  and  finite 
element  predicted  values  was  attributed  to  there  being  only  a  finite 
number  of  crack  extension  paths  available  in  the  finite  element  model 
compared  to  an  infinite  number  in  the  analytical  Si h  theory.  However, 
the  finite  element  model  always  predicted  crack  extension  along  the 
closest  available  crack  extension  path  to  that  of  the  theoretical  crack 
extension  direction. 
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3.2  Off-Axis,  Unidirectional  Composite  Tension  Specimens. 

A  tension  test  of  a  unidirectional  coupon  with  a  small  crack  pre¬ 
sent,  Fig.  (1.3a),  was  simulated  as  a  test  of  the  crack  extension  direc¬ 
tion.  The  mesh  used  consisted  of  306  elements  and  338  node  points.  Fig. 
(D.l).  Referring  to  Fig.  (3.1a)  and  (2.7),  the  boundary  conditions  were 
chosen  to  simulate  the  grips  of  a  tension  test  machine  and  are  given  by: 

at  y  =  0:  v(y  =  0,z)  -  -6,  w(y  =  0,z)  =  0 
at  y  =  2L:  v(y  =  L,z)  =  6,  w(y  =  L,z)  =  0 
at  z  =  0,28:  traction  free 

The  dimensions  were  chosen  to  be  similar  to  a  typical  tensile  coupon, 
(Reference  [19]  suggests  L  >  30B  be  used.),  and  were  taken  as, 

3  =  0.25' 1 
L/2B  =  15 

The  applied  displacement  load  was  chosen  to  be, 

6  ;  0.00!"  =  375o(in)  (3.14c) 

The  initial  half-crack  length  was  chosen  as, 

0.005  in 
sine 


(3.14a) 

(3.14b) 


(3.13a) 

(3.13b) 

(3.13c) 


a 


(3.14d) 
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and  the  other  parameters,  the  fiber  direction,  6,  the  crack  inclination 
angle,  B,  and  the  crack  position.  A,  were  varied  to  simulate  various 
conditions.  The  problem  was  treated  as  being  a  case  of  orthotropic 
plane  stress  with  the  constitutive  relation  being  that  of  Eqn.  (B.13). 
The  material  system  chosen  was  T300/5208  graphite-epoxy,  properties 
listed  in  Appendix  C.  The  crack  was  assumed  to  be  a  virtual  crack, 
i.e.,  having  zero  initial  width. 

3.2.1  Comparison  of  the  Crack-Fxtension  Direction  Theories  for  a  30° 
Lamina 

The  theories  for  predicting  the  crack  extension  direction,  i.e., 
the  modified  Griffith  theory.  Section  (2.2.1),  the  Sih  strain  energy 
density  theory.  Section  (2.2.2),  the  Tsai-Wu  theory.  Section  (2.2.3), 
the  modified  point  stress  theory.  Section  (2.2.4),  and  the  modified 
Hashln  theory.  Section  (2.2.4),  were  compared  against  the  expected  crack 
extension  path  for  a  30°,  off-axis,  unidirectional  tensile  specimen  with 
a  crack  orientated  along  the  z  axis,  i.e.,  6  =  90°  in  Fig.  (3.1a). 

i  . 

The  predicted  crack  extension  path  for  the  modified  Griffith  theory 
is  shown  in  Fig.  (3.3a),  the  Sih  strain  energy  density  criterion  in  Fig. 
(3.3b),  the  Tsai-Wu  theory  in  Fig.  (3.3c),  the  modified  point  stress 
theory  in  Fig.  (3.3d)  and  the  modified  Hashin  theory  in  Fig.  (3.3e). 

The  experimentally  observed  direction  is  shown  in  Fig.  (3.3f).  Experi- 

\ 

mental  results  for  graphite-epo*y  [6,  7  and  8],  Boron-aluminum  [20],  and 
graphite  polyimide  [21],  unidirectional  composites  all  Indicated  that 
failure  of  unnotched  spec-mens  and  fracture  of  notched  specimens 
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occurred  along  planes  parallel  to  the  fibers.  Fig.  (3.3a)  thru  (3.3e) 
indicate  that  the  modified  point  stress  and  the  modified  Hashin  theories 
were  the  only  theories  that  predicted  the  correct  crack  extension  path. 

The  reason  that  the  modified  Griffith,  the  strain  energy  density 
and  the  Tsai-Wu  theories  predicted  incorrect  crack  extension  paths  was 
because  none  of  these  theories  account  for  the  differences  in  the  energy 
required  to  create  crack  extension  surfaces  at  arbitrary  angles  to  the 
fibers. 

Since  the  modified  point  stress  and  Hashin  criteria  were  the  only 
criteria  to  yield  accurate  results,  further  case  studies  were  limited  to 
these  two  criteria. 


3.3.2  Variation  of  the  Modified  Point  Stress  and  Hashin  Functions  for  a 
30°  Lamina. 

The  modified  point  stress  function,  R ( rQ , «|> ) ,  was  given  by  Eqn. 
(2.20)  as 


R(r0,4>) 


(3.15) 


where  was  taken  as  the  normal  stress  and  T  defined  by  Eqn. 

(2.22),  the  Strength  normal  to  a  plane  of  possible  crack  extension.  The 
modified  Hashin  function  can  be  defined  as  H(rQ,4>),  where  from  Eqn. 
(2.33), 


H(r  ,♦)  =  (^)2  +  £*±) 
T.  .  T  . 


(3.16) 
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where  o  was  the  shear  stress  and  T  the  shear  strength,  Eqn.  (2.34), 

r<j>  *  t 

on  the  plane  of  possible  crack  extension. 

Finite  element  predictions  of  R (rQ ,4> )  and  were  computed  as 

a  function  of  4>  for  three  different  values  of  rQ  for  the  30°  lamina  of 

Section  3.2.1.  (The  finite  element  mesh,  the  boundary  conditions  and 

the  crack  geometry  were  identical  to  those  of  the  problem  considered  in 

Section  3.2.1).  The  values  for  the  R(ro»$)  function  were  normalized  to 

their  maximum  values  and  plotted  in  Fig.  (3.4)  for  three  different 

normalized  values  of  r  .  The  rQ  values  were  normalized  with  respect  to 

rmax  where  rmax  was  the  limiting  value  of  r  for  the  finite  element 
o  ’  o  ° 

mesh  used,  (e.g.,  see  Section  2.2.4.4).  The  values  for  the  H(rQ,*) 
function  of  the  modified  Hashin  theory  were  also  normalized  to  their 
maximum  value  and  were  plotted  in  Fig.  (3.5)  for  the  same  three  norma¬ 
lized  r  values. 

0 

The  results  of  Fig.  (3.4)  and  (3.5)  achieve  maximum  values  at  +  = 
210°  for  all  of  the  various  r0/r™ax  values.  Since^the  criteria  both 
assume  crack  extension  in  the  direction  where  R(ro,4>)  or  H(rQ,4>)  reach  a 
maximum,  both  criteria  choose  the  expected  value  of  4>  =  210°,  (e.g.,  see 
Section  3.2.1).  Figure  (3.4)  and  (3.5)  also  indicate  that  the  trends 
remain  unchanged  though  the  values  differ  slightly  for  the  three  rQ/r0 
values.  This  indicates  that  the  prediction  of  the  crack  extension 
direction  is  fairly  insensitive  to  the  value  chosen  for  rQ.  Note  that 
in  the  case  of  the  modified  point  stress  theory.  Fig.  (3.4),  the  maximum 
value  of  R/Rmax  increased  with  increasing  distance  from  the  crack  tip. 
This  trend  was  just  the  opposite  of  what  was  expected  since  the  stresses 


0.0^ 


270 

<J>  (degrees) 


Figure  3.4  Normalized  Modified  Point  Stress  Ratio  vs .<t> 
For  A  30?  Unidirectional,  Laminate 
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are  actually  singular  at  the  crack  tip,  i.e.,  at  rQ  =  0.  The  reason  for 
this  discrepency  is  that  as  rQ/r™ax  approaches  unity,  the  stresses  are 
computed  nearer  the  element  corners  which  introduces  numerical  errors  in 
the  stress  computation.  This  stands  as  a  good  example  of  why  such 
criteria  are  limited  to  the  prediction  of  crack  extension  direction  and 
not  the  load  which  would  cause  crack  extension.  The  values  for  the 
normalized  modified  point  stress  ratio,  Fig.  (3.4)  and  the  normalized 
modified  Hashin  ratio.  Fig.  (3.5),  for  0°  <  *  <  =  180°  were  not  shown 

because  of  numerical  difficulties  in  the  computation  of  the  stresses  in 
this  regime.  However,  the  modified  point  stress  results  seemed  to 
indicate  a  slight  increase  in  R(rQ,4»)  for  0°  <  *  <  30°  and  some  negative 
values  between  *  =  30°  and  *  =  180°.  Similarly  for  the  modified  Hashin 
ratio.  Fig.  (3.5),  the  values  for0°  <  *  <  180°  seemed  to  indicate  slight 
Increases  in  H(rQ,*)  around  *  =  30°  and  *  =  120°  with  some  negative 
values  for  between  *  =  30°  and  +  =  180°.  Comparison  of  the  modified 
point  stress  ratio,  Rfr^),  Eqn.  (3.15),  and  the  modified  Hashin 
ratio,  H(rQ,ij>),  Eqn.  (3.16),  reveal  that  j 

I 

t 

»(V*)  =  CR(r0,O3Z  +  (^L)2  (3.17) 

r«f> 

Hence,  the  differences  in  Fig.  (3.4)  and  (3.5)  represent  the  effects  of 

squaring  the  point  stress  ratio,  R(rQ,*),  plus  the  shear  effect  of 

including  the  square  of  the  shear  ratio,  o  /T 

r<T  r<)>*  l 
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3.2.3  Crack  Extension  in  Unidirectional  Laminates 

The  modified  point  stress  and  the  modified  Hashin  criteria  were 
used  to  predict  the  crack  extension  direction  and  the  modified  crack 
closure  technique  was  used  to  predict  the  energy  release  rate  for 
several  different  unidirectional  laminates.  The  reason  the  study  of 
crack  extension  in  unidirectional  laminates  was  limited  to  the  modified 
point  stress  and  Kashin  criteria  was  because  they  were  the  only  criteria 
which  accounted  for  differences  in  the  energy  required  to  create  new 
fracture  surfaces  at  arbitrary  angles  to  the  material  principal  system-- 
as  demonstrated  for  a  30°  lamina  in  section  3.2.1.  The  mesh,  boundary 
conditions  and  material  system  were  identical  to  those  explained  at  the 
beginning  of  section  3.2.  A  virtual  crack,  i.e.,  one  having  zero  ini¬ 
tial  width,  was  assumed.  The  orientation  of  the  crack  was  assumed  to 
be  B  =  90 •  and  A  =  L/2,  Fig.  (3.1a). 

The  results  are  shown  in  Table  (3.3)  where  0  is  the  material  prin¬ 
ciple  coordinate  system  orientation,  $o  the  orientation  of  the  plane  of 
crack  extension,  G  is  the  total  energy  release  rate,  Eqn.  (2.4),  and  the 
%  mode-I  and  %  mode-II  values  represent  the  mode-I  and  mode-II  contri¬ 
butions  to  the  total  energy  release  rate,  G.  The  expected  crack  exten¬ 
sion  directions  were  based  on  experimental  results  for  graphite-epoxy 
[6,  7],  Boron-aluminum  [20],  and  graphite  polyimide  [21]  composites. 

The  results  of  these  experiments  all  indicated  that  failure  of  unnotched 
specimens  and  fracture  of  notched  specimens  occurred  along  planes  par¬ 
allel  to  the  fibers  and  that  was  always  greater  than  180°  except  in 
the  case  of  6  =  0°,  where  extension  could  occur  in  either  a  $  =0° 


Table  3.3  Finite  Element  Crack  Extension  Predictions 
For  Unidirectional  Laminates 
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or  180°  direction.  The  results  of  Table  (3.3)  indicate  that  both 
the  modified  point  stress  and  Hashin  criteria  choose  the  correct  path  of 
crack  extension  in  all  of  the  cases  considered.  The  fact  that  both 
criteria  also  choose  dual  values  of  <j>o  =  0°,  180°  for  the  0  =  0°  speci¬ 
mens  was  also  promising  since  this  is  what  would  be  expected.  While  no 
direct  comparison  of  the  computed  energy  release  rates  could  be  made, 
the  computed  values  were  similar  to  and  in  the  same  range  as  experimen¬ 
tal  values  obtained  by  Wang  and  Crossman  £63  for  double  side  notched 
graphite-epoxy  specimens.  The  %  mode-I  and  the  %  mode-II  values  were 
also  compared  with  the  results  of  Wang  and  Crossman  [6]  and  were  found 
to  exhibit  the  same  trends  and  range  of  results. 

3,2.4  Effects  of  Crack  Orientation  on  the  Fracture  Characteristics  of 
Unidirectional  Laminates 

The  finite  element  model  was  used  to  predict  the  fracture  charac¬ 
teristics  of  a  30  ,  uni di recti onal  laminate.  The  cases  considered  werer 

j 

j 

1.  Influence  of  crack  position,  referring  to  Fig.  (3.1a),  a  crack  near 

the  grip  of  the  tensile  machine  was  simulated  by  specifying 

that  £  =  0.1  with  B  =  90°. 

L  \ 

2.  Influence  of  specimen  aspect  ratio,  referring  tp  Fig.  (3.1a),  a 

\ 

short  specimen  was  modeled  having  length  L/2B  =  2  with  3  =  90°. 

3.  Influence  of  crack  orientation,  referring  to  Fig!  (3.1a),  specimens 
with  3  =  30°,  3  =  50°,  and  3  =  90°  were  considered. 
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All  of  the  cases  assumed  that  a  virtual  crack  existed  and  used  the  same 
mesh,  boundary  conditions,  loading,  material  and  geometry,  unless  stated 
otherwise  above,  as  described  at  the  beginning  of  section  3.2. 

The  results  of  case  1,  Table  (3.4),  indicated  that  the  crack  direc¬ 
tion  for  the  crack  near  the  grip,  A/I.  =  0.1,  extended  in  the  same  direc¬ 
tion,  4>  ,  as  did  the  center  cracked  specimen,  A/L  =  0.5,  while  the 
o 

energy  release  rates,  6,  and  the  %  mode-I  and  .%  mode-II  contributions  to 
the  total  energy  release  rate  differed  by  10-20%. 

The  case  2  results.  Table  (3.5),  indicated  that  the  crack  extension 
direction,  ,  was  insensitive  to  the  specimen  aspect  ratio  and  that  the 
resulting  energy  release  rates  differed  by  about  20%.  However,  the  type 
of  fracture  that  occurred  was  the  opposite  for  the  two  cases.  That  is, 
the  long  specimen,  L/28  =  15,  fractured  in  a  mainly  mode-I  fashion  while 
the  short  specimen,  L/2B  =  2,  fractured  in  a  mainly  mode-tl  fashion. 

This  trend  is  not  surprising  since  Nemeth,  et  al.  [19]  have  shown  that  a 
completely  different  stress  state  exists  in  short  specimens,  i.e.,  L/2B 
=  5,  as  compared  to  long  specimens,  i.e.,  L/2B  =  15. 

The  results  of  case  3,  Table  (3.6)  predicted  the  same  crack  exten¬ 
sion  directions,  <t>Q,  and  similar  energy  release  rates  and  fracture  modes 
for  the  three  crack  inclination  angles,  3,  considered.  These  results 
are  consistent  with  the  experimental  results  of  [6],  [7],  [20],  and  [21] 
which  indicate  that  fracture  of  unidirectional  laminates  always  occurs 
on  planes  parallel  to  the  fibers. 


Table  3.4  Influence  Of  Crack  Position  On  The  Fracture 

Characteristics  Of  30?  Unidirectional  .  Laminates 


Table  3.5  InfluenceOf  Specimen  Aspect  Ratio  On  The 
Fracture  Of  3( P,  Unidirectional,  Laminates 
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3.3  Free  Edge  Crack  Growth  in  Laminated  Composite  Tensile  Specimens 

Cracks  located  at  the  free  edge  of  laminated  tension  specimens. 

Fig.  (2.10),  were  considered.  The  crack  extension  paths  and  the  result¬ 
ing  energy  release  rate  vs.  crack  length  curves  were  considered. 

The  solution  procedure  consisted  of  obtaining  a  generalized  plane 
strain  solution  of  the  front  face,  i.e.,  the  y-z  plane,  under  an  applied 
load  of  then,  using  the  results  of  the  front  face  solution  as 
applied  loads,  obtain  a  fracture  mechanics  solution  to  the  free  edge, 
i.e.,  the  y'-z*  plane,  via  subsequent  generalized  plane  strain  solu¬ 
tions.  (The  generalized  plane  strain  formulation  is  described  in  Appen¬ 
dix  A.) 

The  generalized  plane  strain  solution  of  the  front  face  was 
obtained  using  a  132  element  by  150  node  finite  element  mesh.  Fig. 

(D.3).  Referring  to  Fig.  (2.10)  and  Eqn.  (A. 4),  the  boundary  conditions 
and  geometry  used  were: 


B  = 

=  n.25" 

H  = 

0.02",  e  =0.001 

XX 

(3.18a) 

at 

y  =  o. 

u(y  = 

0,z)  =  V(y  =  0,z)  -  0  ! 

(3.18b) 

at 

y  -  b, 

tract i 

Ion  free  \  j  . 

(3.18c) 

at 

z  =  0, 

W(y,z 

=  o)  =  o  V 

\ 

(3.18d) 

at 

2  =  H, 

tracti 

on  free 

1 

(3.18e) 

where  quarter  symmetry  was  assumed. 

Since  the  generalized  plane  strain  solution  assumes  stresses  and 
strains  to  be  independent  of  the  out  of  plane  coordinate,  the  study  of 


the  free  edge  crack  growth  problem  can  be  modeled  exactly  only  when 
there  is  no  variation  of  these  quantities  in  both  the  x  and  x'  direc¬ 
tions.  Since  the  x-direction  is  the  out  of  plane  direction.  Fig. 

(3.6a),  for  the  front  face  problem  and  the  x'-direction  is  the  out  of 
plane  direction,  Fig.  (3.6b),  for  the  free  edge  problem,  the  study  is 
limited  to  cases  which  exhibit  x  and  x*  independence.  The  only  lami¬ 
nates  which  exhibit  this  quality  are  unidirectional  laminates  because  of 
the  absence  of  edge  effects.  However,  if  the  thru  the  thickness  varia- 


bles  of  the  free  edge  problem.  Fig.  (3.6b),  are  assumed  constant 

and 

>  ;  w 

i  j  \ 

equal  to  the  values  at  the  free  edge  of  the  front  face  problem,  i 

.e.,  at 

y  =  B  in  Fig.  (3.6a),  then  the  corresponding  free  edge  boundary  condi- 

.  -  -  _ 

tions  become  Fig.  (2.10),  (3.6a),  (3.6b), 

v.  • 

at  y'  =  0,  W (y '  =  0,z')  =  W(y  =  B,z) 

(3.19a) 

V*  (y*  =  0,z' )  =  -[U(y  =  B,z)  +  ^-(x  =  C)] 

(3.19b) 

•  ">  - 

!  k; 

at  y'  =  b,  M'(y'  =  2C,z')  =  W(y  =  B.z) 

(3.19c) 

i  i  • 

V(y*  =  2C,z')  =  -[U(y  =  B.z)  +  exx*(x  =  -C)] 

(3.19d) 

\  \  ’ 

i  ■ 

i 

1  ■; 

j 

where  2C  is  the  free  edge  modeled  length. 

1  \ 

at  z'  =  0,  W (y*  ,z'  =  0)  =  0 

(3.19e) 

:  'N 

at  z'  =  H,  traction  free 

(3. 19f ) 
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with: 

U * (y ' , z ' )  =  0  (3.19g) 

exx  =  ex^z')  =  TT^y  =  B,z) 

and: 

e!  =  6.  -  90° 

Several  laminate  configurations  were  considered  as  candidates  for  the 
analysis  of  free  edge  crack  growth.  Tests  were  preformed  to  find  which 
laminate  configurations  could  best  be  modeled  as  generalized  plane 
strain  problems  on  both  the  front  face  and  the  free  edge.  The  tests 
consisted  of  first  obtaining  the  front  face  solution  then  applying  the 
corresponding  boundary  conditions,  Eqns.  (3,19a-i),  to  an  uncracked  free 
edge  model.  The  results  of  the  tests  were  analyzed  to  determine  which 
laminates  gave  the  best  correspondence  between  the  stresses  and  strains 
at  the  free  edge  as  well  as  approximated  independence  of  the  displace¬ 
ments  in  the  x  and  x'  directions.  Fig.  (2. 1C),  for  the  front  face  and 
free  edge  problems.  The  results  indicated  that  both  angle-ply  and 
cross-ply  laminate  configurations  gave  a  reasonable  correspondence  of 
free  edge  stresses  and  strains.  However,  only  the  cross-ply  config¬ 
urations  approximated  independence  of  both  x  and  x\  Fig.  (2.10),  under 

\ 

the  applied  load.  Hence,  the  analysis  of  free  edge  crack  growth  in  this 

study  was  limited  to  crcss-ply  laminates.  i 

! 

Experimental  results  £6,  22]  indicate  that  laminates  which  contain 
90°  plys  along  with  other  plys  where  0  90°  can  exhibit  transverse 


(3.19h) 
(3.19i ) 
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crack  growth  in  the  90°  plys  and  del ami  nations  between  the  90°  and  0  ^ 
90°  plys.  As  an  example  of  the  quasi-three  dimensional  capabilities  of 
the  finite  element  model  the  nature  of  crack  growth  in  a  [902/0;>j|s, 
T300/5208,  graphite-epoxy  tensile  specimen  was  considered.  Referring  to 
Fig.  (3.7),  both  transverse  crack  growth  and  delamination  crack  growth 
were  considered  and  the  results  of  the  two  analysis  were  compared.  The 

method  of  analysis  consisted  of  testing  two  cases.  Case  1  was  to  trace 

the  crack  growth  of  an  initial  transverse  crack  and  case  2  was  to  trace 

the  crack  growth  of  an  initial  delamination.  Both  cases  used  the  boun¬ 

dary  conditions  of  Eqns.  (3.19a-c).  The  transverse  crack  case  used  the 

j 

crack  geometry  of  Fig.  (3.1b)  and  the  delamination  case  used  the  crack 
geometry  of  Fig.  (3.8).  Both  cases  were  modeled  as  the  free  edge  of  an 
8-ply  tensile  specimen  with 

B  =  0.25in,  H  =  O.Olin,  C  =  2H,  a  =  0.77H  (3.20) 

and  an  applied  normal  strain  of.  Fig.  (2.10), 

exx  =  0.001  (3.21) 

The  modified  point  stress  criterion.  Section  2. 2. 4.1,  was  used  to  pre¬ 
dict  the  direction  of  crack  extension  and  the  modified  crack  closure 
method  was  used  to  compute  the  energy  release  rates. 

The  resulting  crack  growth  sequence  for  the  initial  transverse 
crack,  case  1,  is  presented  in  Fig.  (3.9).  The  case  1  results  indicated 
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that  the  transverse  crack  extended  through  the  90°  plys  down  to  the  0° 
plys  when  it  turned  into  a  delamination.  Note  that  this  does  not  neces¬ 
sarily  indicate  that  a  transverse  crack  will  turn  into  a  delamination 
but  just  that,  if  the  transverse  crack  was  the  only  crack  present  and  if 
it  was  to  continue  to  extend  that  it  would  extend  as  a  delamination 
before  extending  through  the  6=0°  plys.  The  plot  of  the  corresponding 
energy  release  rate  vs.  crack  length  is  shown  in  Fig.  (3.10).  Point  A 
represents  the  initial  crack  length  and  point  B  represents  the  point  at 
which  the  transverse  crack  turned  into  a  delamination.  The  critical 
energy  release  rate  for  such  a  crack,  is  given  by  Wang  and  Crossman  [6] 
for  a  mode-I  transverse  crack  in  graphite-epoxy  as 

GIC  -  O.^-1!5-  (3.22) 

in£ 

Since  the  available  energy  release  rate  curve.  Fig.  (3.10),  was  less 

than  0.6  ^n,g-s-  for  the  entire  crack  extension  sequence,  the  crack  would 
in 

not  have  extended  at  the  applied  load  level  of  =  0.001.  However,  if 
the  load  were  increased  the  G  vs.  a  curve  would  have  transletea  up, 
retaining  the  exact  same  shape, [6],  C-D  in  Fig.  (3.10)  for  example,  and 
would  have  predicted  crack  growth  until  the  crack  length  was  such  that 
the  available  energy  release  fell  below  the  critical  energy  release 
rate.  Regardless  of  applied  load,  however,  the  shape  of  the  G  vs.  a 
curve.  Fig.  (3.10),  indicates  that  G  is  decreasing  with  increasing  a. 
Hence,  there  exists  a  certain  range  of  e  values  for  which  the  trans¬ 
verse  crack  would  grow  then  arrest  when  G  falls  below  Gc.  This  trend  is 


G  (in  lbs 


Figure  3.10  G  vs.  a  For  Transverse  Crack  Growth 
In  A  Cross-Ply  Laminate 


supported  by  many  experimental  results  of  [6,  2?],  for  example,  which 
show  that  transverse  cracks  will  grow  through  the  90°  plys,  arrest  and 
form  other  transverse  cracks  at  regular  spacings. 

The  results  for  the  initial  delamination  crack,  case  2,  indicated 
that  the  initial  del  ami  nation  would  immediately  turn  into  a  transverse 
mode/  In  comparing  the  point  stress  ratios  of  the  initial  delamination 
crack  with  the  initial  transverse  crack  at  the  point  it  turned  into  a 
delamination,  it  was  found  that  the  point  stress  ratio  was  roughly  two 
times  larger  for  the  delamination  turning  into  a  transverse  crack  than 
for  the  transverse  crack  turning  into  a  delamination.  This  indicates 
that  if  both  a  delamination  and  a  transverse  crack  were  present  in  a 
given  specimen,  the  delamination  would  turn  and  grow  into  a  transverse 
crack  before  the  initial  transverse  crack  would  extendo  Th>$  trend  is 
supported  by  the  experimental  results  of  Harris  and  Orringer  [22],  which 
indicate  that  transverse  cracks  can  branch  off  from  delaminations. 


Chapter  4 
CONCLUSIONS 

The  present  investigation  has  been  concerned  with  predicting  the 
direction  of  crack  extension  as  well  as  the  load  to  cause  extension  in 
composite  materials.  The  results  of  the  present  study  indicate  that 
failure  criteria  can  be  used  to  predict  the  direction  of  crack  extension 
and  that  an  energy  release  rate  approach,  implemented  through  the  use  of 
a  modified  crack  closure  integral,  can  be  used  to  determine  when  a  crack 
extends  and  if  crack  arrest  will  occur.  The  finite  element  models 
presented  herein  were  formulated  for  two-dimensional  and  quasi  three- 
dimensional  analysis.  However,  the  procedures  and  methods  developed  can 
be  applied  to  full  three-dimensional  analyses  as  well. 

It  was  found  that  criteria  for  predicting  the  direction  of  crack 
extension  should  account  for  differences  in  the  energy  required  to 
create  crack  surfaces  at  arbitrary  angles  to  the  material  principal 
system.  The  Griffith  criterion,  the  Tsai-Wu  failure  criterion  and  the 
Sih  strain  energy  density  theory  all  were  unsatisfactory  in  this  regard; 
but  the  modified  point  stress  and  Hashin  criteria  provided  good  predic¬ 
tions  for  crack  growth  direction.  It  was  also  found  that  since  the 
procedure  developed  herein  assumes  that  crack  extension  will  occur  along 
the  element  sides  adjacent  to  the  crack  tip  node,  an  incorrect  direction 
of  extension  can  be  chosen  if  no  element  sides  coincide  with  the  actual 
direction  of  extension.  However,  the  model  always  chooses  the  closest 
direction  available  to  the  actual  direction  of  extension.  This  problem 


88 


6'  -1 


i 


89 

can  easily  be  overcome  by  first  finding  an  approximate  direction  of 
extension  using  a  coarse  mesh  then  refining  tho  mesh  in  the  area  of  the 
direction  chosen  by  the  coarse  mesh. 

The  present  investigation  has  also  shown  that  future  research  is 
warranted  in  the  following  areas: 

1.  The  extension  of  the  current  work  into  a  full  three-dimensional 
model.  A  full  three-dimensional  model  should  be  formulated  to 
account  for  the  three-dimensional  crack  growth  characteristics  of 
many  laminates. 

2.  Experimental  work  in  the  area  of  critical  e~  vgy  release  rates. 
Experimental  work  should  be  performed  to  determine  the  effect  on  the 
critical  energy  release  rate  of  cracks  extending  at  arbitrary  angles 
to  the  material  principal  system  so  that  the  load  to  cause  failure 
can  be  accurately  defined. 
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APPENDIX  A 


Linear  Elastic  Finite  Element  Relationships 

The  finite  element  used  is  shown  in  Fig.  (A.l).  The  element  is  a 
four  node,  general  quadralateral ,  isoparametric  element.  The  element 
uses  linear  interpolation  as  described  by  Segerlind  [23].  The  details 
of  the  finite  element  concept  are  also  given  in  reference  [23],  The 
technique  involves  mapping  a  distorted  shape  in  the  Cartesian  (y,z) 
Coordinate  System  into  a  square  in  the  Local  (£,n)  Coordinate  System 
where  £  and  n  range  from  ~1  to  +1.  The  relationship  between  the  global 
Cartesian  and  the  local  coordinates  is 

Y  =  N1Y1  +  N2Y2  +  N3Y3  +  N4Y4  =  NiYi  i  =  1,4 

Z  =  NjZx  +  N2Z2  +  N3Z3  +  N4Z4  =  NiZ1  1  -  1,4  (A.l) 

where  the  N. (£,  n)  are  the  interpolation  functions  for  the  four  node 
points  and  Y^,Z;  are  the  Cartesian  coordinates  of  the  nodes.  The  inter¬ 
polation  functions  are  given  by  ! 

I  [  ' 

Nj  -  }(l-€)(l-n),  N2  =  |(l+C)(l-nj 

N3  =  |(lH)(l+n),  N4  =  i(l-0(l+n)  \  (A.2) 

\ 

For  an  isoparametric  element,  the  same  interpolation  functions  are  used 

| 

for  the  assumed  displacements  as  for  the  geometry.  Hence,  for  plane 

93  I 


95 


stress  or  plane  strain,  whe-e  there  exists  two  unknown  degrees  of  free¬ 
dom  per  node,  the  interpolation  becomes 


v(y,z)  =  H1vi  i  =  1,4 

w(y,z)  -  NiWi  1-1,4  (A. 3) 

For  generalized  plane  strain,  where  there  are  three  degrees  of  freedom 
per  node,  the  interpolation  is  given  by 

u(x,y,z)  =  U(y,z)  +  exx.x  -  +  e^.x  1  =  1,4 

v(x,y,z)  =  V(y,z)  =  N<»1  1  =  1,4  (A.4) 

w(x,y,z)  =  W(y,z)  =  j  =  1>4 

where  u,  v  and  w  are  the  x,  y  and  z  displacements,  respectively,  u^ ,  v^ 
and  wf  are  the  unknown  values  at  the  ith  node  and  are  functions  of  y  and 
z  only,  and  is  the  total  strain  in  the  x  direction,  which  is  assumed 
to  be  constant  and  either  known  or  unknown. 

The  strain-displacement  relationships  are  derived  based  on  the 
small  strain  -  small  displacement  theory.  For  the  three  dimensional 
(generalized  plane  strain)  case,  these  relationships  may  be  written  as 


96 


ORIGINAL  PAGc  t** 

OF  POOR  QUALITY 


f 

— 

— 

r 

e 

3 

0 

0 

XX 

3x 

G 

0 

3 

0 

yy 

3y 

r  m 

£ 

ZZ 

0 

0 

3 

3"z 

u 

‘V 

0 

3 

3z 

3 

3y 

V  * 

\z 

3 

3z 

0 

3 

3x 

w 

Yxy 

3 

3 

3x 

0 

^  J 

L 

— 

(A.5) 


And  for  plane  strain  or  plane  stress  these  relationships  become 


yy 


zz 


yz 


a 

3y 

0 

3 

3z 


0 

3 

az 

3 

3y 


Substituting  (A. 4)  into  (A.5),  for  generalized  plane  strain. 


(A.6) 


e 

XX 

h  0  0 

e 

yy 

0  ^y  0 

e 

ZZ 

0  0  h 

o  2-  2_ 

y 

yz 

u  3z  3y 

3  n  3 

y*z 

JI  °  Tx 

3  3  n 

Yxy 

_3y  3x  °_ 

CIN1|IN2|IN3|IN4]{q} 


A 

where  , 

I  is  the  3x3  identity  matrix  j 

fq}  is  the  12  x  1.  vector  of  nodal  displacements 
gi ven  by 


xx 

0 

0 

0 

0 

0 


(A.7) 


.  'V 
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{q} 


U1 

V1 

W1 

u2 

v2 

w2 

> 

Uo 

v3 

w3 

u4 
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(A.8) 


Substituting  (A. 3)  into  (A. 6),  for  plane  stress  or  plane  strain. 


where 


r  *1 

C 

3 

0 

o 

yy 

3y 

e 

zz 

0 

3 

az 

>> 

— ..  * 

3 

TI 

3 

W 

1 

- 

ClN1!lN2|lN3|IN4]{q} 


I  is  the  2  x 
{q }  is  the  8  x 


2  identity  matrix 

1  vector  of  nodal  displacements  given  by 


(A.9) 
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{q}  = 


•vf 

wi 

v2 

w2 

v3 

w3 

v4 

_w4. 


(A. 10) 


The  [B]  matrix  (strain-displacement  relationships)  for  generalized  plane 
strain  can  be  defined  as 


{e}  =  CB]{q] 

(6x1)  (6x12) (12x1) 


(A. 11) 


where  by  comparing  (A. 11)  and  (A. 7),  for  generalized  plane  strain. 
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[B]  = 


0 

0 

3 

0 

^y 

0 

3 

7z 

3 

3 

Jz 

W 

0 

3 

3 

3x 

0 

[IN1|IN2|IN3|IN4] 


(A. 12) 


For  plane  strain  or  plane  stress,  the  [B]  matrix  can  be  writtei.  as 


M  =  CB]{q} 
(3x1)  (3x8) (8x1) 


(A. 13) 


Comparing  (A. 13)  and  (A. 9),  for  plane  stress  or  plane  strain. 


[B]  = 


[IN1|IN2|IN3|IN4] 


(A. 14) 


Recall  from  (A. 2)  that  the  are  functions  of  the  local  coordinates 
£  and  n.  In  order  to  determine  the  elements  of  the  [B]  matrix  a  rela¬ 
tionship  between  the  derivatives  in  the  global  (y,z)  and  the  local 
(£,  n)  coordinate  systems  is  needed.  This  relationship  is  given  by  the 
Jacobian  matrix  [J]  of  the  transformation  where; 
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and 


[0]  = 


L3n 


3N.n 

'3N.‘ 

*  -  CO]-1  ■ 

3£ 

3N. 

1 

3Ni 

Tz 

w 

3z1 

7T 


3z 

3n 
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(i  -  1,  2,  3,  4) 


(A. 15a) 


(A. 15b) 


Substituting  (A.l)  into  (A. 15) 


[J] 


8N1  9N2  3N3  3N4 

W  TC  W  W 


3N,  3N„  3N-  3N 


1  '"’2 
TiT  TrT 


_ 3  __4 

3n  3n 


Yl 

Y2 


Z1 

Zo 


*3  z3 
y4  z4 


(A. 16) 


The  stress-strain  relationships  are  derived  in  Appendix  B.  For 
general  purposes  they  can  be  written  as 


{o)  =  [0]{c} 


(A. 17) 


For  generalized  plane  strain,  { a }  and  {c}  are  6x1  vectors,  and  are 
given  in  Appendix  B  by  Eqns.  (B.9)  and  (B.10),  respectively. 

For  generalized  plane  strain,  the  [D]  matrix  takes  on  the  values  of 
the  [C]  matrix  as  given  by  Eqn.  (B.8).  For  plane  stress  or  plane 
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strain  {a}  and  { r. }  are  the  3  x  1  vectors  given  by  Eqns.  (B.14)  and 
(B.15),  respectively.  For  orthotropic  plane  stress  the  [D]  matrix  takes 
on  the  value  of  the  3x3  [0]  matrix  of  Eqn.  (B.13),  while  for  isotropic 
plane  strain  the  [D]  matrix  becomes  the  3  x  3  [C]  matrix  of  Eqn.  (B.18). 

The  total  potential  energy,  *,  of  a  given  finite  element  Is  the  sum 
of  the  strain  energy,  U,  and  the  work  of  external  loads,  W.  The  strain 


energy,  U,  of  the  element  is 


I)  -  |  J  Mte|Td  * 

V 

and  the  work  of  external  loads,  W,  is 
u  =  -(q)|f}T 


(A. 18) 


(A. 19) 


where  {F J  is  appl ied  mechanical  load,  (traction),  vector. 

Hence,  the  total  potential  energy  of  the  element  is  given  by 


1.  =  U  +  W 


(A. 20) 


Substitution  into  (A. 20)  In  terms  of  the  matricies  and  vectors  described 

\  \ 

herein  and  minimizing  with  respect  to  the  unknowns  yields  the  finite 

\ 

element  equations 


KHq}  -  {Fl 


(A.21) 


\ 


i 
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where  [K]  is  the  elemental  stiffness  matrix  given  by 

[K]  =  J  [B]T[D][B]dV  (A.22) 

V 

when  niapped  into  the  coordinate  systero 

dV  =  |J|d£dn  (A.23) 

where  |J|  is  the  determinate  of  the  Jacobian  matrix.  The  limits  of 
integration  are  -1  to  +1  in  both  5  and  n.  Hence,  Eqn.  (A.22)  becomes 

[K]  =  J  J  [B]T[D][B]|J|dCdn  =  /  /  G(5,n)d5dn  (A. 24) 

-1  -1  -1  -1 

In  order  to  evaluate  the  stiffness  matrix,  [K],  a  numerical  integration 
is  necessary.  Using  a  2  x  2  Gauss  rule  Eqn.  (A. 24)  can  be  evaluated  as 

2  2 

[K]  =  Z  Z  H  H  G(a  ,b  ) 
i=l  j=l  1  J  J 

where  G(£,n)  =  [B]T[D][B]| J | 

(a-,b)  are  the  coordinates  of  the  four  Gauss  points  given  by 
^  J 


(A. 25) 

(A. 26) 
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/ 


and  are  the  corresponding  weight  functions  given  by 


^  =  1,1 

Hj  =  1.1 


(A. 28) 


The  elemental  matricies  are  then  assembled  into  a  global  system  of 
equations,  the  prescribed  boundary  conditions  are  imposed  and  the  system 
of  equations  is  solved  for  the  unknown  displacements.  (This  procedure 
is  explained  in  most  finite  element  books  [23].) 

The  strains  can  be  found  at  any  C,n  location  within  an  element 
through  the  use  of  Eqn.  (A. 7)  for  generalized  plane  strain  and  through 
Eqn.  (A. 9)  for  plane  stress  or  plane  strain.  The  strains  can  be  conver¬ 
ted  to  stresses  by  using  the  stress  strain  relation  of  Eqn.  (A. 17)  where 
[D]  takes  on  the  values  of  the  appropriate  constitutive  relation  (Appen¬ 
dix  B). 
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APPENDIX  B 

CONSTITUTIVE  RELATIONS 

The  const' tutive  relation  for  an  orthotropic  material  in  the  prin 
cipal  material  directions.  Fig.  (B.l),  is  given  by  Jones  [24}  as 


{°)l  «  CC]{e}j 


(B.l) 


where 


[C]  = 


C11  C12  C 


13 


C22  C23 
c33 


(symmetric) 


0 

0 

0 


C44  0 

c, 


0 
0 
0 
0 

55  0 
C66 


(B.2) 


Wi- 


'11 


22 


33 


23 

13 


12 
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A 


(B.3) 
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(B.4) 


1-v  V 

23  32 


‘  ’  Cl2  ’  EUE33S  ’  ‘  E11E22‘ 


V  +V  V 

12  32  13 


.  cn 


V  +V  V 

13  12  23 


1-V  V 
1  12  21 


r  1-v13v31  ,  V23w2l'’l3  c  . 

23  "  EnE33S  ’  23  "  E11E22‘  33 


(B.5) 


lil 


C44  s  G23»  C55  =  613»  C66  =  612 


A  = 


1'V12V21~V23V32~V31V13~2V21V32V13 

E11E22E33 


and 


^31 

E33 


(B.6) 
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For  a  6  rotation  about  the  3,  (z),  axis.  Fig.  (8.1),  the  consti¬ 
tutive  relation  becomes 


{*} 

■  mw 

where 

Cn  c12 

^13 

0 

0 

C16 

c22 

C23 

0 

0 

C26 

[C]  » 

C33 

0 

0 

C36 

(Symmetric)  C44  £45  0 

C65  0 


(B.10) 
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If  m  =  cos  0  and  n  =  sin  0  then  the  components  of  the  [C]  matrix  are 
given  by 


4  2  2  4 

C 1 1  =  m  Cjj  +  2m  n  (Cj2+2Cgg)  +  n  ^22 

=  m  n  ^  +  ^  +n  ^12 

2  2 

C13  =  rn  C13  +  n  C23 

C16  =  mn[m  (Cn-Ci2-2C66)  +  n  (C12~C22+2C66^ 
^22  ~  ^  ^11  +  ^  ( 0 1 2"^ 2C g 6 ^  ^22 

^23  =  n  C13  +  m  C23 

C26  =  mn^-n  ^CH"C12_2C66^+ni  ^C12-C22+2C66^ 

C33  =  C33 

C36  =  mn(C13-C23) 

2  2 

^44  =  m  C44  +  n  C55 
C45  =  mn(C55-C44) 

^55  =  n?C44  +  m  C55 

^66  =  m  n  (C j \~^\2+^"2Z^  ^66^  ^ 


(B.ll) 


For  an  orthotropic  material  under  a  state  of  plane  stress  with 
a  0  rotation  about  the  3,  (x)  axis.  Fig.  (B.2),  Jones  [24]  gives  the 


constitutive  relation  as 


1°}  =  CQ3{e} 


(B.12) 
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Figure  B.2  Material  And  Global  Coordinate  System 
For  A  Two-Dimensional  Lamina 


I 

I 
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APPENDIX  C 


MATERIAL  PROPERTIES 

The  material  properties  for  graphite-epoxy  T300/5208  are  given  by  Nagar- 
kar  and  Herakovich  [25]  as 


Elastic  Moduli 

En  =  19.2  x  106  psi,  E22  =  1.56  x  106  pSi,  E33  =  1.36  x  10^  pSi 

Shear  Moduli 

G23  =  0.487  x  106  psi,  G13  =  0.820  x  106  psi,  G12  =  0.820  x  106  psi 

Possion's  Ratios 

v23  =  0.490,  v13  =  0.238,  v  =  0.238 

6.35  x  103  psi 
23.85  x  103  psi 
12.6  x  103  psi 


Strength  Parameters 

Xj  =  219.5  x  103  psi,  Yy  =  6.35  x  103  psi,  Zy  = 

Xc  =  -246.0  x  103  psi,  Yc  =  -23.85  x  103  psi,  Zc  =  - 

S23  =  9.8  x  103  psi,  S13  =  12.6  x  103  psi,  Sj2  = 


The  material  properties  for  the  isotropic  problems  considered  were 
chosen  to  be 

E  =  30  x  106  psi,  v  =  0.20,  =  50  x  103  psi 

where  ou  is  the  ultimate  strength  of  the  material. 
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APPENDIX  D 


FINITE  ELEMEi'T  MESHES 
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APPENDIX  E 


ON  THE  RELATIONSHIP  BETWEEN  FREE-EDGE  STRESSES  AND 
THE  DIRECTION  OF  CRACK  EXTENSION  IN  ANGLE-PLY  LAMINATES 

It  was  originally  thought  that  the  tree-edge  stress  state  could  be 
used  to  predetermine  the  direction  of  crack  extension. 

Herakovich  [7,  8],  for  example,  has  demonstrated  that  the  crack 
types  present  on  the  free-edge  of  tension  specimens  seem  to  exhibit 
distinct  crack  growth  patterns  direction  of  extension  for  different 
angle-ply  laminates.  To  account  for  this  effect  it  was  thought  that  the 
orientation  of  the  principal  stress  plane  on  the  free  edge  would  follow 
along  the  same  path  as  the  crack  plane  when  traced  through  the  laminate 
thickness.  To  test  this  premise,  the  orientation  of  the  crack 
plane,  4>  ,  through  the  laminate  thickness  was  compared  with  the  through 

H 

the  thickness  variation  of  the  principal  stress  plane.  A  generalized 
plane  strain  solution  of  the  front  face,  Fig.  (3.6a),  was  obtained  using 
the  same  geometry,  loading  and  boundary  conditions  as  given  in  Section 
3.3.  Only  the  stresses  in  the  plane  of  the  free  edge  were  considered, 
i*e*>  °xx*  °zz  and  t  ,  in  Fig.  (2. IP).  It  was  assumed  that  the  thermal 
curing  stresses  were  balanced  by  the  hygroscopic  stresses  in  the  lami¬ 
nate  so  that  only  mechanical  loading  be  considered.  The  material  con¬ 
sidered  and  experimental  crack  plane  orientation  was  that  of  Herakovich 
[7]  (the  material  was  the  same  T300/520S  graphite-epo>y  of  this  report— 
with  properties  listed  In  Appendix  C).  The  orientation  of  the  principal 
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plane,  <J>  ,  was  calculated  through 

r 

♦  =  270°  -  lun-'Cjfjf-;]  (E.l) 

XX  zz 

The  results  for  four  laminate  configurations  are  presented  in  Fig. 
(E.la-d).  The  multi-valued  angles,  4>p,  at  z/H  =  0.5  represent  the 
effect  of  a  crack  turning  into  a  delamination  at  that  point.  Fig. 
(E.la-d)  indicate  that  only  in  the  [902/02]s  case  did  the  theory  agree 
with  the  experiment.  Hence,  it  was  concluded  that  the  direction  of 
crack  extension  could  not  be  predicted  through  the  crack  free  stress 
state  but  that  an  actual  crack  must  be  introduced. 


\ 
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APPENDIX  F 


COMPOSITE  LAMINATE  FRACTURE  MECHANICS  (CLFM2D) 
INPUT  DATA  SEQUENCE 


Cards  1 

&  2:  title  cards  (20A4) 

Column 

Contents 

1-BO  ITITLE  (I,J) 

Title 

Card  3: 

Control  card 

(11IS) 

Column 

Contents 

1-5 

NPROB 

Problem  type  (1  =  generalized  plane 

6-10 

NEM 

strain,  2  =  skewed  plane  stress, 

3  =  orthotropic  plane  stress,  4  =  iso¬ 
tropic  plane  strain) 

Number  of  elements  in  mesh 

11-15 

NODS 

Number  of  nodes  in  mesh 

16-20 

NANG 

Number  of  different  angles  r  must  be 

>  1 

21-25 

NSDF 

Number  of  specified  degrees  of  freedom 

26-30 

NSBF 

Number  of  specified  forces  (tractions) 

31-35 

NEXX. 

Number  of  different  Exx  -  must  be  >  1. 
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36-40  NFX  Fracture  mechanics  key 

NFX  =  0,  no  fracture  analysis 
NFX  f  0,  fracture  mechanics  problem 
41-45  NPLOT  Plot  key 

j  NPLOT  =  0,  no  plot 

NPLOT  =  1,  undeformed  plot  only 
NPLOT  =  2,  deformed  plot  only 
NPLOT  =  3,  both  deformed  and 
undeformed  plots 

46-50  NT1  Oump  file  for  displacements 

NT1  =  0,  no  dump 
NT1  f  0,  dumps  displacements  to 
unit  NT1 

51-55  NT2  Dump  file  for  strains,  stresses  and 

strain  energies 

NT2  =  0,  no  dump 

NT2  i  0,  dumps  to  unit  NT2 

56-60  NCHECK  Data  check  option  (NCHECK  f  0  -  data 

i 

check  only) 

\ 


Card  4: 

Scale  factors 

(2D10.5) 

Column 

Contents 

1-10 

SCAY 

Y-scale  factor 

9 

11-20 

SCAZ 

Z-scale  factor 
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Card  5: 

Printer  control 

card  (SI5) 

Column 

Contents 

1-5 

KEY  (1) 

Key  for  printing  element  data 

6-10 

KEY  (2) 

Key  for  printing  nodal  data 

11-15 

KEY  (3) 

Key  for  printing  specified  displacement 

data 

16-20 

KEY  (4) 

Key  for  printing  specified  force  data 

21-25 

KEY  (5) 

Key  for  printing  nodal  displacements 

26-30 

KEY  (6) 

Key  for  printing  strains,  stresses  and 

strain  energies 

Note:  If  KEY  (I)  f  0,  print 

Card  6:  Plotter  scale  factors  (4D10.5)  -  skip  if  NPLOT  =  0 
Column  Contents 


1-10 

PVSCL 

Plot  Y-scale  factor 

11-20 

PZSCL 

Plot  Z-scale  factor 

21-30 

VMAX 

Maximum  Y-displacement 

31-40 

WMAX 

Maximum  Z-displacement 

Card  7: 

Material 

property  card  1  (6D10.5) 

Column 

Contents 

1-10 

PROP  (1) 

E11 

11-20 

PROP  (2) 

e22 

21-30 

PROP  (3) 

e33 

31-40 

PROP  (4) 

G23 
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41-50 

PROP  (5) 

g13 

51-60 

PROP  (6) 

G12 

Card  8: 

Material  property  card  2  (6D10.5) 

Column 

Contents 

1-10 

PROP  (7) 

v23 

11-20 

PROP  (8) 

v13 

21-30 

PROP  (9) 

°12 

31-40 

PROP  (10) 

“11 

41-50 

PROP  (11) 

“22 

51-60 

PROP  (11) 

“33 

Card  9: 

Angles  (8D10.5) 

Column 

Contents 

1-10 

ANG  (1) 

Angle  No.  1  (in  degrees) 

11-20 

• 

ANG  (2) 

Angle  No.  2  (in  degrees) 

• 

• 

• 

ANG  (NANG) 

• 

Angle  No.  "NANG"  (in  degrees) 

124 


Card  10 

:  Element  data  (5x, 

1415)  repeat  “NEM"  times 

Column 

Contents 

1-5 

Blank  (can  put  in  element  numbers) 

6-10 

NOD  (1,1) 

Node  No.  1  of  element  I* 

11-15 

NOD  (1,2) 

Node  No.  2  of  element  I1 

16-20 

NOD  (1,3) 

Node  No.  3  of  element  I* 

21-25 

NOD  (1,4) 

Node  No.  4  of  element  I* 

26-30 

IANG(I) 

Angle  number  of  element  1^ 

31-35 

IEPS  (1,1) 

O 

exx  number  at  local  node  1 

36-40 

IEPS  (1,2) 

3 

exx  number  at  local  node  2 

41-45 

IEPS  (1,3) 

O 

exx  number  at  local  node  3J 

46-50 

IEPS  (1,4) 

3 

exx  number  at  local  node  4 

51-55 

ISTRS  (1,1) 

Stress  output  location  l4 

56-60 

ISTRS  (1,2) 

Stress  output  location  24 

61-65 

ISTRS  (1,3) 

Stress  output  location  34 

66-70 

ISTRS  (1,4) 

Stress  output  location  44 

t  . 

71-75 

ISTRS  (1,5) 

Stress  output  location  5^ 

l 


1.  Refer  to  Fig.  (F.l)  for  node  numbering  sequence. 

2.  Refer  to  Fig.  (F.2)  for  angle  orientation. 

3.  exx  numbers  correspond  to  those  of  Card  16. 

4.  Stress  output  locations  shown  in  Fig.  (F.l). 

| 

i 

s 

i 

i 

i 

i 

| 

I 


i 

\ 


127 


ORIGINAL  PACL  »3 
OF  POOR  QUALITY 


Card  11:  Nodal  data  (5x,2D10.5)  -  repeat  "NODS"  times. 


Column 

Contents 

1-5 

Blank  (can  put  in  node  numbers) 

6-15 

Y(I) 

Y  coordinate  of  node  I 

16-25 

Z(I) 

Z  coordinate  of  node  I 

Card  12:  Specified  displacement  data  (115,010.5)  -  repeat  "NSDF" 
and  skip  if  NSDF  =  0. 

Column  Contents 


1-5 

ND 

Node  number 

6-10 

IDR 

Direction  (1  =  x,  2  =  y,  3  = 

11-20 

UBDF(I) 

Specified  displacement  value 

Card  13 

:  Specified  force 

data  (215,  DIO. 5)  -  repeat  "NSBF 

and  skip  if 

NSBF 

=  o. 

Column 

1 

Contents ! 

1-5 

ND 

Node  number 

6-10 

IDR 

Direction  (1  =  x,  2  =  y,  and 

\ 

11-20 

UVSF(I) 

Specified  force  value 

\ 

Card  15: 

:  Temperature 

data 

■  ■■  '  V\ 

(1010.5)  1 

Column 

Contents  j 

1-10 

TEMP 

Temperature  change 

times 
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Card  16:  Normal  (x-di rection)  strain  data  (8D10.5) 
Column  Contents 


1-10 

EPSX(l) 

eyx  value  No.  1 

11-20 

EPSX(2) 

e  value  No.  2 

A  A 

• 

# 

# 

• 

» 

• 

EPSX(NEXX) 

exx  value  No.  "NEXX" 

*1  f  NFX 

=  0,  end  of  data 

Card  17:  Fracture  control 

card  (615,2010.5) 

Column 

Contents 

1-5 

NODE 

Initial  crack  tip  node  number* 

6-10 

N0DE1 

Secondary  node  No.  1* 

11-15 

N0DE2 

Secondary  node  No.  2* 

16-20 

NSE 

Number  of  stop  elements  (must  be  >  1) 

21-25  - 

NSN 

Number  of  skip  nodes  (must  be  >  1) 

26-30 

NCRT 

Fracture  criteria 

NCRT  =  1  Griffith 

NCRT  =  2  Sih-strain  energy  density 

NCRT  =  3  3-0  point  stress 

NCRT  =  4  2-D  point  stress 

31-40 

ANOT 

Initial  crack  length 

41-50 

TMAX 

Maximum  run  time  (CPU-seconds) 

1.  Refer  to  Fig.  (F.3)  for  crack  nodes  description. 


ORJGtJyAL  WA&ti; 

OF  POOR  QUALITY 


NODE  1  NODE2 


Figure  F.3  Nodas  Defining  Crack  Tip 
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* 


( 


Card  18:  Stop  element 
Column 

1-5  NSTOP  (1) 

6-10  NSTOP  (2) 


data1  (1615) 

Contents 

Stop  element  No.  x 
Stop  element  No.  2 


;  \ 


I 

| 

i 

( 


NSTOP  (NSE)  Stop  element  No.  "NSE" 


Card  19:  Skip  node  data^ 
Column 

1-5  NSKIP  (1) 

6-10  NSKIP  (2) 


(1615) 

Contents 
Skip  No.  1 
Skip  node  No.  2 


NSkIP  (NSN)  Skip  node  No.  "NSN" 

1.  Stop  elements  terminate  the  solution  when  they  are  reached  - 
used  to  prevent  tear  through. 

2.  Skip  nodes  are  used  to  eliminate  the  designated  nodes  as 
crack  growth  possibilities  -  used  to  prevent  tear  through 
and  the  crack  from  growing  back  on  itself. 


\ 


; 

5 


4 
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Card  20 

:  Strength  properties  (3010. 5) 

Column 

Contents 

1-10 

STRENG  (1) 

*T 

11-20 

STRENG  (2) 

yt 

21-30 

STRENG  (3) 

h 

end  of  data 


